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ANALYSIS OF SURFACE ABLATION OF NONCHARRING MATERIALS 


WITH DESCRIPTION OF ASSOCIATED COMPUTING PROGRAM 

By Fred W. Matting and Dean R. Chapman 
Ames Research Center 


SUMMARY 


A generalized method is presented for solving the problem of stagnation- 
point heat transfer and material response for blunt bodies experiencing melt- 
ing and vaporizing or subliming ablation. The analysis is applicable to 
wind-tunnel and flight conditions (with body forces taken into account) ; inter- 
nal radiation can be considered or the material can be assumed opaque; the 
analysis can be used for different planets. During entry flights, a body will 
start in the free-molecule regime, pass through a transitional regime, and 
finish in the continuum regime of gas dynamics. Approximate equations, 
rationally obtained, are presented which provide "bridges" between the free- 
molecule and continuum regimes for such quantities as the convective heating 
rate, surface shear, heat blockage, and mass loss. Several illustrative exam- 
ples (including surface chemical reaction cases) show the applications of the 
analysis. Comparisons with experiment are made where possible. The analysis 
has been machine programmed for numerical solutions using a finite difference 
scheme, and a running energy balance is kept as a check on accuracy. Instruc- 
tions are provided for using the computing program. 


INTRODUCTION 


The development of heat shields for space vehicles and long-range 
missiles has motivated a marked increase in the effort to understand the pro- 
cess of ablation. This intensified study is contributing to the understanding 
of natural ablative phenomena that occur when extraterrestrial bodies enter 
the Earth* s atmosphere. Ablation data obtained in the laboratory, for example, 
in arc -jet wind tunnels, do not duplicate in any single experiment all the 
conditions of entry flights; hence, there is need for analytical methods of 
predicting and explaining the ablative phenomena. Before such methods can be 
applied, however, their validity and accuracy must be determined by comparing 
calculated results with results from wind-tunnel tests, with flight data, or 
with post-flight observations of a man-made or natural object when one can be 
recovered. 

A generalized method is presented here for solving the combined problem 
of heat transfer and material response for the stagnation region of blunt 
bodies experiencing melting and vaporizing or subliming ablation. An attempt 
has been made to describe the problem mathematically as completely as possible 
in order to obtain nearly exact solutions. This required that the analysis be 
machine programmed for numerical solutions. This program in various stages of 



development has heen used successfully at Ames Research Center during the 
past three years. It has previously heen employed in the analysis of tektite 
ablation (refs. 1, 2). 

A number of the equations used in the analysis are already well 
established (refs. 3-6); however, some of the equations and features of the 
analysis are new. The analysis takes account of the fact that entry bodies 
initially fly in the free -molecule regime, then in a transitional regime, and 
finally in the continuum regime of gas dynamics. The analysis contains for- 
mulas for the transitional regime to bridge between the free -molecule and 
continuum regimes; these formulas have been rationally derived from simple 
models and are believed to fill an important gap in previous analyses of small 
objects entering a planetary atmosphere. 

Several options in the analysis and associated computing program are 
available. Internal radiation in the body is accounted for, or the body can 
be assumed to be opaque. Flight cases as well as wind-tunnel cases can be 
calculated; the flight cases can be applied to any planet, provided certain 
characteristics of the atmosphere are known. The rear boundary conditions for 
the ablating material can be those for a heat shield, or the aerodynamic base 
heating for an object, such as a tektite, can be accounted for. The ablating 
material can be a type that melts and/or vaporizes, sublimes, or undergoes a 
surface chemical reaction in the ablation process. The various material prop- 
erties and the external flow conditions can be put into the computing program 
arbitrarily so that a variety of ablation research problems can be studied. 

To illustrate the types of problems that can be handled by the analysis 
(and to elucidate some of the significant ablative phenomena), several typical 
examples are presented. These are calculated by the numerical computing pro- 
gram associated with the analysis. The principal features of the numerical 
computing program are described, and instructions are given for its use. 


ANALYSIS AND METHOD OF SOLUTION 


In this section, the principal emphasis is on the method of solution 
used in the associated computer program. Equations are presented in specific 
form with numerical values inserted for universal constants. In performing 
calculations for various materials and situations, the user has a number of 
parameters (constants) at his disposal; typical values of these are listed in 
the description as they appear. 


Basic Approach and Approximations 

The analysis is concerned with the problem of surface ablation near the 
stagnation point for transparent or opaque materials of finite thickness. The 
analysis is thus restricted to materials that undergo surface ablation includ- 
ing melting, evaporation, sublimation, or surface chemical reactions, such as 
depolymerization. Chemical reactions involving the external gases are not 
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considered. The "boundary-layer equations, as such, are not solved, but 
results of known solutions are used to obtain heating rates, surface shear, 
and the effects of mass transfer. Front and back surface heating rates must 
be solved for, in order to determine the front and back face boundary condi- 
tions. Two types of heating environments are treated. One is of constant 
velocity and stream density, such as generally exists in an arc-heated wind 
tunnel; the other is of time-varying velocity and density over an entry flight 
(using a constant lift/drag ratio but a variable mass). In this latter envi- 
ronment, the equations of motion are solved simultaneously with the heat- 
transfer and ablation equations. In the equations of motion, the ratio of 
body mass to the product of drag coefficient by frontal area appears, and this 
is related (empirically) to the surface recession. 


Conservation Equations 


The analysis is essentially a time -dependent energy balance along the 
stagnation center line of an axisymmetric blunt body. The basic equations to 
be solved are the conservation equations for energy, mass, and momentum for 
the material of the body, written in a simplified form that is valid in the 
stagnation region. The energy and momentum equations have been further sim- 
plified by neglecting inertial terms; this procedure is a valid approximation 
for viscous ablative materials, such as glass, stone, or any subliming mate- 
rial. The curvilinear coordinate system used is shown in sketch (a). The 

conservation equations are : 

SHOCK 



Sketch (a) 


Energy 


/ST _ ST 

PC Ijt + V ^ 


_S_ 

Sy 


K ^ - F(y,t) 


( 1 ) 


where F(y, t) is the internal radiation 
flux term. 


Continuity (for constant density) 

u + Sv = 0 

Sx x dy 

x Momentum 

3 f 5u\ _ Sp _ pax 
3y v dyj “ Bx R 


( 2 ) 


( 3 ) 


Most previous analyses have taken account of internal radiation or vehi- 
cle acceleration, but not both. The acceleration term in the momentum equa- 
tion (3) can be important in determining the rate of removal of melted 
material; this term, as written, is valid in the stagnation region. 
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The continuity and momentum equations can be put into quadrature form. 

We make use of the fact that _u varies linearly with x near the stagnation 
center line, and we define: u = (Su/5x) x=0 . Then, near the center line 

u = ux ( 4 ) 

The continuity equation in quadrature form is: 


u dy 1 (5) 

The value of v w depends on surface temperature and external conditions; its 
evaluation is described below in the section. Front Face Velocity. The sur- 
face recession rate and surface recession are, respectively. 



v sr = H = |VBF| = K^bF'^I 



(6a) 

(6b) 


Pressure Gradient 

The solution to equation (5) can be substituted into equation (l), but to 
solve equation ( 5 ) we must first obtain u(y,t) from the momentum equation. 

We write, for the pressure near the center line. 


P = p(0) + | p"(o)x 2 


(7) 


and equation ( 3 ) becomes, after substituting equation ( 4 ), 

f ) - p,,(o) - ? ■ - p " <8) 

where P ,T is then the negative of the second derivative of pressure with a 
correction for the body force acting on the material. We can integrate equa- 
tion ( 8 ) twice and obtain: 


u(y,t) = P" 


>y- 


J y 


BF y 1 dy 1 + , 

| J. w 


f: 


BF dy^ 
M- 


(9) 


where t -^- 1 (dynes/cm 3 ) is the x derivative of the absolute value of the 
surface shear, (dT^/dx) x==0 . In p" the quantity, a, is the acceleration of a 
body in flight (positive for increasing speed) . For wind-tunnel calculations 


4 



a = -g for a vertical wind tunnel with, upward flow and a = 0 for a 
horizontal wind tunnel. In the calculations to he made, P " (dynes/cm 4 ) is 
evaluated as follows: 

Wind tunnel: 


Flight: 


«•* = (a. - M 


P" = 


PR 


. i M \ 

W. 


4mr 2 


2 x 10 4 DV c 
R 2 


K. 


tu 


where A 2 (approximately unity) is a correction used with a modified 
Newtonian approximation. 


(10a) 


(10b) 


= Voo [ 2A^ 

to R 7 P 21 


(11) 


and Eis in equation (lOa) is selected to give the wind-tunnel body force. 
The quantity, K-^ u , is an empirical correction for oscillation or tumbling in 
flight. It is evaluated as 

K t u = K + (1 - K) 


1 - e 


-(X/Xi) 


E 


16 


( 12 ) 


where K is the fractional time that the center point is initially exposed to 
approximately stagnation conditions, and and Ei 6 are values selected to 
give a realistic damping to the oscillation during entry. Typical values for 
the constants that have been used for calculating a tumbling tektite entry are: 
K = 0 . 25 ; ^*1 = R /10 cm; Ei 6 = 1 . 0 . For a nontumbling body, K is unity and 
the other constants have no influence . 


Surface Convective Heat Transfer 


In order to evaluate the x gradient of wall shear, t v * , we require 
first an evaluation of the surface convective heat transfer; we also need this 
to determine the surface boundary conditions. Instead of enthalpy (cal/g), 
we use an ''enthalpy velocity" (km/sec), which is defined as 

V 2 = O.OO836 h s (13a) 


V = 0.0915 -/h^ 


(13b) 
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For the external gas we use an average specific heat. 


c 


P 


1 

T 



( 14 ) 


In laminar continuum flow we can evaluate the surface convective heat transfer 
as (ref. 7 ) (with a vorticity correction) 


< 3 oc = A 4 71 vl,15 (v 2 - O.OO836 c p T w ) (1 + (15) 

where A 4 is a constant and C 6 is a vorticity correction (generally small). 
The results of equation (15) agree well with existing experimental data over 
an extended range of enthalpy potentials. The value of A 4 for Earth entries 
is approximately 1 . 1 ; in wind-tunnel tests, A 4 is evaluated with a calorim- 
eter. With a blowing correction we have 

V = H>c ( l6 ) 

where \|r is evaluated in equation (28) . 

For high altitude flight or for rarefied wind-tunnel conditions, it is 
necessary to consider the free-molecule regime. For surface convection in 
the free-molecule regime we use a Newtonian type of approximation (ref. 8, 

pp. 395-403). 


4 FM - ^ ■ 0 - 00 ® 3 4 * 6 ( 17 ) 

In the transition regime between free-molecule and continuum flow, the con- 

vective heat transfer will have a value bridged between the evaluations in 
equations (l 6 ) and (17)- This has been derived from a simple kinetic theory 
model (see appendix C) yielding the result 


%w 



( 18 ) 


When the tumbling correction is applied, the convective heat transfer at the 
front wall is 


^w ^w^tu (^ 9 ) 

Equations (15) to (19) are used to determine the convective heat transfer 
to the front face, but it is also of interest to calculate some related quan- 
tities. When no material is being lost to the vapor state, \|r = 1 , and 
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equation (l8) specializes to 


^oo 



e 




( 20 ) 


and equation ( 19 ) specializes to 




( 21 ) 


When material is lost to the vapor state , ^ ^ 1; we can consider that \|r 
operates on q Q c but not on qpM* Then a modified \|/ can be obtained to 
operate on Q . We define the modified blowing factor, \jF, as: 

f = 

^OO 

^■\|fW ^o (22b) 

On multiplying both sides of equation (22b) by K^, we have also 



a w = ^ 


(23) 


When__we substitute equations (l 6), (l8), and (20) into (22a), the evaluation 
of \|; becomes 

t A _ e -W*<loc\ 
t = 

1 - e 



Equations (20) to (24) are alternate forms entirely equivalent to equa- 
tions (l8) and (l9)* Tbe quantities q^, q Q , and \]F, although not needed in 
determining the convective heat transfer, are computed as quantities of 
interest conceptually. 

The bridging relation given in equation (l8) or the alternate forms in 
equations (20) and (24) will automatically take account of changes of heating 
rates as a body flies from one regime into another and will place control in 
the appropriate regime. Comparisons with available measured data are shown 
in appendix C. 

For normal (nonrarefied) wind-tunnel conditions, the bridging relations 
given above are not needed; also there is no tumbling, so we have simply: 
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Ktu = 1 


(25a) 




(25b) 


*oc 

(25c) 



(25a) 


In the heat-transfer relations given above from equation ( 15 ) to equa- 
tion (25), the blowing parameter, is needed. To evaluate t , we need to 
know the relationship between the equilibrium vapor pressure and temperature 
for the ablating material, p ve = p ve (T w ). We will write 


P = 


Pt- 


- 1 = 


• vm 






L*ve 




- 1 


(26) 


where p m is the modified equilibrium vapor pressure when the equilibrium is 
shifted by the presence of other gaseous materials in the boundary layer. An 
example of this is the suppression of vaporization of silica due to the pres- 
ence of oxygen in the atmosphere. This suppression effect for silica is 
analyzed in reference 9 in which an analytic, but implicit, expression is 
obtained for the modified equilibrium vapor pressure. The use of the exponent, 
E 7, is an empirical accounting for this effect which yields values within 
several percent of the values obtained by the rigorous analysis of reference 9* 
For silica, the value E 7 = 1.4 has been used. For a number of other mate- 
rials, E 7 = 1 , or the modified equilibrium vapor pressure is the usual equi- 
librium value. For p-t 2 (in atmospheres) we use a hypersonic approximation 
(twice the free-stream dynamic pressure with a correction factor): 


AJV 2 

P t2 101-3 


( 27 ) 


where a typical value for Ai is 0.95j and the number, 101. 3, accounts for 
the units in the equation. Equation (26) is a limiting form for continuum 
conditions; P could be based on the actual existing vapor pressure, p v , 
which will not be an equilibrium value, but in the limit p v approaches p m 
from below. However, p v cannot exceed p^ , so f° r ‘this equation, whatever 
value p-^ may actually take, it is not allowed to exceed pt 2 in equa- 
tion (26). This has been done by giving P a lower limit which is a small 
positive number (lO -6 ). So P is defined by equation ( 2 6) down to its lower 
limit. The quantities, p v , p ve , p w , p^ 2 , and P, are considered to be 
evaluated at the liquid or solid surface. 

We can represent for an evaporation or sublimation process as 

f - -n 35 + E35 (28a) 

1 + V 

p 



This relation, with the asymptote, E 35 = 0 . 06 , gives a good fit to a number of 
boundary -layer solutions (refs. 10 - 12 ). Equation ( 28 a) is also a limiting 
form for continuum conditions since P, as evaluated in equation (26), is not 
based on the actual vapor pressure, p v . The quantity l/P is a function of 
the mass loss rate due to vaporization (see eqs . (38), ( 3 9a), and ( 40 ) below). 
For a surface chemical reaction we use the form 


t = 


1 - E fl5 

^ B11 

+ P [1 + ( v wc/ v wFM) ] 


+ S35 


(28b) 


where v WG and are si 1011 ™- evaluated velow in equations (38) and (39e). 

The correction using v wc and v w fm takes account of the possible reaction 
rate control of the mass loss rate which does not depend on the (modified) 
equilibrium vapor pressure (eqs. ( 39 e) and (kO ) ) . The quantity, B u in equa- 
tion (28) (0. 95-1. 55 used), depends on the ratio of molecular weights of exter- 
nal gas to blowing gas and should preferably be determined by experiment. In 
the absence of experiment, we can estimate Bu as follows. We can define 



m v 


and we have 


where the constant 


Bn 


constant 

— n 


M 


0.7 to 0.8 and n ^ 2/3 to 3/^* 


(29) 


( 30 ) 


Wall Shear Gradient 

We are now in position to evaluate the x gradient of the wall shear, 
t ¥ ! = (dT w /dx) x -o« We first consider the case of continuum flow with no blow- 
ing to evaluate t q ! = (dT 0 /dx) x=0 . Using a modified form of Reynolds 
analogy, we can write: 


q-oc 

Ah 


Karo’ 



(31) 


where K a is a constant that depends on the Prandtl number and is unity when 
the Prandtl number is unity. Using equation (ll) we obtain 


A _ ^-2 Tp T R (Ah) 

3 K a 10 5 <loc V oo 


(32a) 
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v = (32b) 

R(&h)\lp2± 

■where A 3 « 1.45 is a typical value. We apply the tumbling or oscillation 
correction, K*t u , (eq. (12)) to t q ? and obtain 


oc 


= V K tu 


(33) 


For the effect of blowing on wall shear in continuum flow we use (ref. ll) 




(34) 


The value of E s should be small; unless experimental evidence is available, 
it is suggested that the value zero be used. (Reasonable answers have also 
been obtained with E 8 = 0.3 B u .) In the numerical computing program, the 
value of P that is inserted into equation (34) (only) is arbitrarily pre- 
vented from becoming less than E 8 so that the shear blowing factor cannot 
become unrealistically large when the actual P is very small. Using equa- 
tions ( 32 b, 33 , and 34), we obtain our expression for t wc * in continuum 
f low. 


T wc * “ 


8 3 6 A 3 Po c VooKtu^ [ ^ - + - X®.©/? ) y ] 
(V 2 - O.OO 836 c p T w ) 


(35) 


For the x 
by blowing ) , we 


gradient of shear in free molecule flow (which is unaffected 
have (ref. 8 , pp. 395-403) 


t _ A cm DVoo 2 10 4 K tu 

t wFM " R 


(36) 


where A cm is the x momentum accommodation coefficient. For the bridging 
between free -molecule and continuum wall shear, we use essentially the same 
model as that used for heat-transfer bridging (see appendix C). 


'w 



./t 1 

'wc 


wFM/ wc 


( 37 ) 


The evaluations of P" and t w ’ given above along with knowledge of the 
temperature distribution enable us to solve equation ( 9 ) for TT(y,t) which 
can be substituted into equation ( 5 ). 
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Front Face Velocity 


To solve equation (5) we need to evaluate the front face velocity, v w . 
For the front face velocity under continuum conditions, we use the so-called 
Lewis analogy (Le = l) , which states that the ratio of mass diffusion to con- 
centration "gradient" is equal to the ratio of continuum heat transfer to 
enthalpy potential: 


0,00836 tq oc 
pMP(V 2 - O.OO836 c p T w ) 


( 38 ) 


An equivalent form of equation (38) is given as equations (31) and (32) of 
reference 9 (see also ref. 13 ) . Equation (38) is a limiting form for contin- 
uum conditions, because of the way P is evaluated in equation (26). (As 
noted previously, P is evaluated at the solid or liquid surface.) Equa- 
tion (38) also contains empiricism in the evaluation of \|r (eq. (28)). In 
reality, it is expected that the diffusion rate should reach a maximum value 
with a very small P, if P were based on the actual pressure at the surface, 
p . This would also require that the i|r asymptote he zero. With a finite 
\|r asymptote (which seems to fit existing data), the calculated diffusion 
rate becomes unrealistically large for small P. Under these conditions it 
can be expected that free -molecule or reaction-rate control will generally 
prevail (see eq. ( 40 )), so that an inaccurate calculation of the diffusion 
rate for these conditions will have little effect on the net rate calculated. 


For the calculation of the front face velocity in the free -molecule or 
rate -controlled regime, we distinguish two cases: (l) evaporation or sublima- 

tion, and ( 2 ) a chemical reaction such as a depolymerization. For the evapo- 
ration or sublimation case, we have, from kinetic theory, the Langmuir 
equation (Cl) . With constants evaluated to account for our units we have: 


v wFMl - 


44.3AcyPve 

P J T w 


( 39 a) 


Using equation ( 29 ) we write 


l v wFMl - 239 


me 


^ve 


29- 1 / P 


w 


( 39 b) 


We now use 


and get 


A f 

A cv 


- A cvJ 


me 

29.1 


| v wFm| 


239Aey^Pye 
p n/mTw 


( 39 c) 

( 39 d) 
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Equation (39^) is the form we use in our calculations. When the molecular 
weight ]% of the external gas is 29 *1* the accommodation coefficient can he 
used directly in equation (39d) ; otherwise the accommodation coefficient must 
he corrected according to equation (39c). 

For the chemical reaction case^ we use an Arrhenius rate form for a first 
order reaction 


v wFMl " Be 


~Et/1w 


(39e) 


In the coordinate system used, both vvc and will he negative quantities. 

The bridging equation between the free -molecule or reaction-rate controlled 
regime and the continuum or diffusion controlled regime turns out to be (see 
appendix C) : 


1 - = + 


v w 


VwFM v wc 


(40a) 


- _ v wFM v wc 
Vw ^wFM + ^wc 

The use of equation (40) automatically places the front face velocity in the 
appropriate controlling regime: the diffusion-controlled, the rate -controlled, 

or the transitional regime. 



Internal Radiation 

In the energy equation (l), we require the evaluation of the internal 
radiative flux, F(y, t) . In the evaluation of F, we assume either an opaque 
body (F = 0) or a transparent gray body (two shades of gray; one for the 
incoming gas cap radiation, another for the absorption -emission of the heat- 
shield material) . A third alternative is to approximate the transparent gray 
body by a semitransparent body that is opaque internally, but has a variable 
surface emissivity. 

The evaluation of the radiative flux for the gray transparent body is 
given below; it is similar to that of reference 14 which treats scattered 
radiation and uses an exponential attenuation. The present evaluation con- 
siders one reflection from the front and the rear surfaces, which is a good 
approximation for materials that absorb well. 
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F(y,t) 

2n 2 acr 


■y yBF 

f T 4 (T))E 2 [a(y - r| ) 1 d.T| - / T 4 (Ti)E 2 [a(ri - y)]drj 

J 0 y 

o^BF 

+ R e ff J T 4 (ri){E 2 [a(y + ti) ] - E 2 [ a (2y BF - (y + T)))]} dr\ 


2n 2 acr 


e" a2y - R eff e 




(41) 


R e ff is the effective coefficient of reflection for planar radiation and is 
related to the maxiraam emissivity hy the relationship (ref. 15 ) 


n R _ max. 

1 “ R eff ~ n 2 


(42) 


The second-degree exponential integral E 2 (z), used as the attenuation func- 
tion in the integrals of equation (4l), is (ref. 1 6 , appendix I): 

/ °° -ZXi pi / 

e Xi 2 ~ dXi = J e Z/|ai d[i 1 ( 43 ) 

In actual numerical calculations, equation (4l) is used directly to evaluate 
F at the front and back faces only. The quantity, g = 8F/8y, is obtained by 
taking the derivative of equation (4l), and g is numerically evaluated for 
use in equation (l) . Internal values of F are then obtained by numerical 
quadrature of the flux derivative, g. The gas cap radiative flux, q^, is 
evaluated with an empirical approximation: 


q R = E 4 PD E5 V E ®EC tu 


(44) 


The form of equation (44) is deduced from experimental correlations presented 
in reference 17- Input constants E 4 , E 5 , and E 6 for a given environment 
should be selected to fit available data. The level of radiation and the 
surface reflectivity are both accounted for in the evaluation of E 4 ; E 5 can 
vary from 0.5 for nonequilibrium radiation to 1.7 for equilibrium radiation 
in air, while E 6 can vary from 5 to 8 for air. After the exponents are 
selected, E 4 should be chosen to give the proper level of radiation. Values 
of the constants that have been used for Earth tektite entries are: 

e 4 = 0 . 76 x 10 “®; e 5 = 0 . 5 ; e s = 7 . 0 . 

In the calculation for a semitransparent body, the material is considered 
to be internally opaque, and the front surface emissivity is varied in an 
appropriate manner with the thermal thickness of the temperature profile. 

This variation is derived by assuming an exponential temperature distribution 
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near the wall. This permits a closed form integration, and the result is 
further approximated to the following simple form: 

€ max _ € max /. ^ 

e FF* ^ 

aA A 

where E 33 = 2.4/a and the thermal thickness. A, is evaluated in equa- 
tion ( 63 ) • (if (dT/dy) v >0 or A < 0, then epp is assigned the value of 
€ max’) With this representation, the hack-face emissivity egp is assumed to 
he constant and F = 0. The semitransparent approximation gives virtually 
the same results for most ablation characteristics as the transparent case 
(see table I), although the internal temperature profiles do not agree closely. 
It greatly reduces computing machine time, however; hence it is used for most 
calculations when an accuracy of the order of 10 percent is adequate. 


Boundary Conditions 

The evaluations of all the terms in the energy equation (l) have been 
shown, so this equation is in a form to be solved for l(y,t) . Boundary con- 
ditions are needed for equation (l), and these are determined in the standard 
manner by writing surface energy balances for the front and back surfaces, 
providing for the appropriate differences between the opaque and transparent 
cases. Options are provided in determining the rear surface boundary 
conditions as shown below. 

The energy balance for the front surface is written: 

- ( K §y) = phv? w + + d - B is)(q. R - 6 FF aT w 4 ) (k6) 


where 

Bie = 0 for the opaque and semitransparent cases 
Bi6 = 1 Tor the transparent case 

In the coordinate system used, the front face remains at the origin. As 
material is lost, the location of the back face recedes to smaller values of 
y^p, or the region in which we are solving equation (l) becomes smaller. So, 
to know the location of the back face, at any time, t, the system of equations 
mast have been solved up to time, t. 

In writing the back surface energy balance, we distinguish between a 
back surface exposed to flight conditions and a back surface in contact with 
a backing material. For the back surface exposed to flight conditions (as, 
for example, with a tektite), we write the following back surface energy 
balance: 
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q BF + (B ls - l) €BF a(T BF 4 - T 0 4 ) 



(47) 


where Biq is zero or unity as noted above. The back surface convective heat 
transfer, q-gp, may be a function of surface temperature. We relate qpp to 
the front surface convective heat transfer (without blowing), q Q , empirically 
as follows: 


q BF _ %/F q o 
(h s - h BF ) (h s - h w ) 


1 + (Ei 9 


l)(l - 


(48) 


where Rp/F is defined as the ratio of base to front -face laminar convective 
heating when normalized by the respective enthalpy potentials. To account for 
transition to turbulent flow at the base, we use: 


%3F turbulent 
^BF laminar 


(49) 


E 


13 - 


0.2 

u Re , . , . 

oo transition 


0.2 


(P10~ 6 ) (Vool0 5 )2Rh 


(J^cA) 


^transition 


1 

at transition 

(50) 


Typical values of constants that have been used for tektite entries are: 

Rr/f = 0.01; E 13 = 0.01; E 19 = 5-0. When the quantities, Rp/p and e-g-p, are 
assigned the value zero, the back boundary condition is adiabatic. 

For numerical computations of the transparent case, the computing program 
has been arranged so that it is possible to modify equations (46) and (47) by 
assigning some radiant energy to the surface energy balances. The reasons for 
this are explained in the section, METHOD OF THE NUMERICAL PROGRAM, Boundary 
Conditions for the Transparent Case, equations (97) an d (98) • 

The other back boundary condition of interest is used for heat -shield 
calculations. We assume a backing material which acts as a heat sink and is 
at uniform (but increasing) temperature, Tpp. For an opaque heat shield the 
heat transfer is entirely by conduction to the backing material; for a trans- 
parent heat shield, we assume that, in addition to conduction, the radiative 
flux, F ( y-g p , t ) is entirely absorbed by the backing material. Then, in place 
of equations (47) to ( 50 ) we write 


c 




+ B 16 F(y^,t) 


(51) 


where c is the heat capacity per unit area of the backing material. In the 
limiting case with c = 0, the back boundary condition is adiabatic, and 
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equation (51) should be replaced by equations (47) and (48) -with %/p and 
e BF assigned the value zero. 

It should be noted that the surface temperatures generally cannot be 
specified a priori, as they depend partly on external conditions that may 
change with time. The computing program must "find" the appropriate boundary 
conditions that satisfy the partial differential equation (l) and the surface 
energy balances. 


Trajectory Equations 

For flight we must solve simultaneously the conservation equations and 
the trajectory equations of motion. We use the two-dimensional trajectory 
equations (with variable mass) for entry in a meridional plane, as shown below 
(ref. 18). 


V oo = V v 2 (52a) 

7 = tan" 1 (52b) 


du 

dt 



uv 

Hr, 


dv _ %> 

dt - R p 10 s 


DV a 


v 


20 


M 


CdA 



dD _ -Dv 
dt Sft 


(53a) 


(53b) 


(54) 


We also use a hypersonic approximation that considers the ambient atmospheric 
enthalpy to have a constant value: 


v = JvJ* + E 3 8 

(55a) 

E 38 = 0.00836 

(55b) 


where h^ is the effective average constant atmospheric enthalpy in cal/g, 
and E 38 (an input to the computing program) has the units km 2 /sec 2 . For 
Earth entries, 0.5 has been used for E 3a . 
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Equations (53) and ( 54 ) can be solved numerically in a straightforward 
manner. As programmed, L/b r has been taken as constant for a given body. 

The atmospheric scale height, S^, has been programmed so that it may change at 
selected values of atmospheric density, D. One computing option for Earth 
entries provides automatic changes in through three values to represent 

the ARDC atmosphere. For arbitrary planet entries, four arbitrary successive 
values of can be used. (See appendix D.) 

The quantity, M/CpA (g/cm 2 ), must be evaluated for use in equation ( 53 )* 
The variation of M/A has been set up empirically as a function of the sur- 
face recession, X. This is equivalent to assuming a geometry for the 
recession shape: 


M 

A 



1 + EiqX 4* C 3 




(56) 


The variation of the drag coefficient, C*q, through the free -molecule, transi- 
tional, and continuum regimes, is represented by a bridging equation developed 
in appendix C: 


C D " C DC 


1 + E 9 e- 15(RD)(l+El4) 


where 

„ C DFM - C DC 

Es ‘ — cS 


(57) 

(58) 


and E 9 depends on the body shape. The free-molecule drag coefficient, 

SpFMj ma y given the value 2 . The parameter, E 14 , has some dependence on 
body shape (and flight conditions), but will often be assigned the value zero 
(the value for a sphere in air). We combine equations (56) and ( 57 ) to obtain 


M _ / M \ 
C D A V C DC A /i 


1 + Ei a X + C3 




1 + Ese -«(BD)(1+K 14 ) 


( 59 ) 


In this equation, (>qc is shown grouped with the initial value of M/Cj^qA, 
because this initial quantity is used in the computing program. The empiri- 
cism in equation (59) can also be used to account for any change in CpQ with 
change of body shape. 


Miscellaneous Relations 

Nose radius .- The effective nose radius, R, is calculated as a quantity 
that has an empirical variation with the surface recession. 
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R = R ± 


(6o) 


1 + e 17 x 2 


+ Ci 




It is desirable to evaluate the constants in equation (60) from experimental 
data. Otherwise, one must estimate the geometrical shape of the ablation 
surfaces. 


Density ratio across normal shock . - In several of the relations presented, 
the density ratio across a normal shock, p 2 i, appears. For wind-tunnel calcu- 
lations, p 2 i is considered to remain constant for a particular case. The 
value of P21 is a portion of the input data for these cases. For flight 
calculations, p 21 changes and must be continuously calculated. For flight 
cases we use the equations: 
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P2 




e -E 12 RV 5 D 



1 + 0.08 log 

°io 




(61) 
(6 2 ) 


Equations (6l) and (62) are empirical relations that give good fits to data 
for air. Equation (6l) accounts for nonequilibrium effects, while equa- 
tion (62) fits equilibrium air data, such as presented in reference 19- The 
equations provide a valid approximation for any gas mixture that consists 
predominantly of nitrogen. The entire analysis is not overly sensitive to the 
evaluation of P 21 * It appears as a square root in equations (ll) and (35) 
which are approximations themselves. Values of En and E i2 that have been 
used for Earth entries are 1.0 and 0.0001, respectively. 


Other Calculated Quantities 

Although not always required in the analysis described in the preceding 
subsections, several other quantities are calculated because they are of 
interest. These are described below. 

Thermal thicknesses . - We calculate the thermal thickness. A, based on the 
wall temperature gradient. 


A = -(Ty - To) 

(8T/8y) w 


(' 63 ) 


This is meaningful only when (dT/dy) w < 0 (otherwise A is arbitrarily 
assigned the value 10 6 ) . We also calculate a viscosity thickness, A^, which 
we define as the depth at which \i = eq^. The corresponding viscosity thick- 
ness temperature is T^. quantity, Ta , is determined from Tw and the 

viscosity representation formula, equation ^85) , below. The quantity, A^, is 
determined by interpolation of temperature profile data. If (8T/8y) w > 0 or if 
Ta^ exists nowhere in the temperature profile, A^ is assigned the value zero. 
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Stored energy comparison with exponential temperature profile . - Using a 
quantity, cp, we compare each calculated temperature profile with an exponen- 
tial profile having the same (ST/Sy) w - The quantity cp is defined as the 
ratio of energy stored to the stored energy associated with the corresponding 
exponential temperature profile (with form similar to equation (Al) in 
appendix A) . Both energies are calculated as constant property approxima- 
tions. For the area under the exponential profile we have: 

Aexp = A(T W - T 0 ) (l - e + R B/F ) „ A(T W - T 0 ) 


Then, 



(6k) 


Removal of melted jmterial by uressuye_gradient and surface shear .- We 
calculate - a” term that compares the removal of melt by the pressure gradient 
and the surface shear. We define the quantity, A 2B , as an approximate evalua- 
tion of twice the ratio of the portion of the melting velocity, 3x, due to the 
pressure gradient to that due to surface shear. The exact quantity would be 
twice the ratio of the first term to the second term in equation (9) > at 
y = 0. However, we use the approximation, p « in the integrals and 

obtain the simple approximate expression, 


A 2D = 


gp'vy 

V 


(6 5 ) 


Surface recession due to vaporiz ation and melting .- We are interested in 
comparing the surface recession due to vaporization with that due to melting. 
We define F as the instantaneous ratio of the surface recession rate due to 
vaporization to the total surface recession rate (due to vaporization and 
melting) . 


F = 



( 66 ) 


For the ratio of the surface recession due to vaporization to the total sur- 
face recession at any time, t: 



Flo w lines in the material . - We are interested in the path of the flow of 
material up to the point of melting or vaporizing. This is of particular 
interest in the study of tektites. To obtain the flow lines, we make an 
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approximate quasi -steady state analysis. 

dx = ux dt 

dy = v dt 
dy/dx = v/ux 


x r (y,t) 


*(y>t) 

X BF 



(u/v)dy 1 


( 68 ) 


where Xgp is any selected small value of x at y B p (where u = 0) . Then 
x = X r xpp gives x(y) at time t for a flow line. 


Aerodynamic deceleration .- For flight cases, we are interested in the 
aerodynamic deceleration of the body. We calculate the ratio of the aero- 
dynamic force to the mass of the body to obtain the aerodynamic deceleration. 
We define apg as the component of deceleration due to aerodynamic drag, 
normalized by the gravitational acceleration of the planet: 


a Dg " 


DV„ 




The absolute value of the total normalized aerodynamic acceleration is 


( 69 ) 


a g = a Dg *7l + ( L/D r ) 2 


( 70 ) 


Reradiation and apparent emissivity .- The terms for the fluxes of 
radiation at the front and back surfaces are given in equation ( 76 ) . The rate 
of reradiation, from the front surface (absolute value) is F^g. For opaque 
bodies. 


f RS - 6 FF°' T w 4 ( 71 ) 

The expression for the reradiation from the back surface of opaque bodies has 
a similar form. For opaque material, the surface emissivities, and e-g-p, 
will be assigned values. For semitransparent material representation, the 
evaluation of e-p-p is as given in equation (45)* 

For transparent bodies, the situation is somewhat different. Here, we 
wish to calculate effective surface emissivities (which will vary with time) . 
The flux of reradiation from the front surface is given by 

Frs = ~ F(0,t) (72) 


Then we calculate: 
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( 73 ) 


_ f rs 

e FF crT w 4 


Similarly, for the "back face: 

F(y BF ,t) 
e BF " oT bf 4 


( 74 ) 


The quantities, and in equations (73) and (7*0 are the calculated 

effective surface emissivities for transparent material. 


Energy Balance 

Equation (l) is solved by numerical finite difference methods , and it is 
desirable to check the accuracy of solutions obtained. This is done by cal- 
culating a group of energy-integral terms listed below, summing them up, and 
determining the residual (error) in the sum. The magnitude of the energy 
integrals is aiso of interest, because it shows the disposition of the 
energies involved; this knowledge gives an insight into the processes of 
ablation. 

We make the listing of the energy rate terms as follows: 

The total convective heat -transfer rate into the material is 

<W = + (T5) 

The net radiative heating rate into the material is 

^-rad “ F(0,t) - ^^BF*^ + ^ B is) - ^FF^^ 4 ~ G BF 0 ^ T BF " -^o 4 ^ 

(76) 

where 

B 16 = 1 for the transparent case 

B 16 = 0 for the opaque and semitransparent cases 

and 

F = 0 for the opaque and semitransparent cases 

The energy accounted for by the rate of vaporization (positive into the mate- 
rial) is 

irap = ph v-v (7T) 
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The energies put into the material are accounted for by terms that involve 
convection of the ablating material in the x and y directions and by a 
storage term. For the rate of increase of stored energy. 


_d_ r y BF 

/^(yj 

Sstor p dt J 

/ C P dT 2 


J T 

L ^o -J 


ay. 

The x direction convection energy rate term (positive out) will be 

>T(y x ) 


r^BF - 

Vcon * 2p V 
u o 


C-p dT2 


T„ 


ay. 


The y direction convection energy rate term (positive out) is 


**w 


Von = c P dTl 

'-'m 

O 

The error in the energy rate balance will be a residual term, l res * 
Sres ^con + ^-rad + ^-vap ^-stor ^-ucon ^vcon 


(78) 


( 79 ) 


(so) 


(81) 


The residual, q res , shows the accuracy of the energy rate balance at any time, 

t • 

The cumulative energy balance shows the total size of the various terms 
involved and the error accumulation. We compute the following integrals: 


Q, 


con 


t i 


<W dtl 


^rad / ^-rad dtl 


Q- 


-vap 


^vap 


^stor p 


y BF( t ) 





pi( y^t) 

pY-op(ti) 

p^(y , t ± ) 

/ C P dT 2 

d y x - p / 

/ <=p d T 2 
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1 

J 0 
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( 82 a) 

( 82 b) 

( 82 c) 

V 

( 82 d) 
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r t 

^acon J ^ucon ^3 


Q. 


vcon 


■£ 


^vcon ^1 


( 82 e) 

( 82 f) 


So, for the accumulated residual we have: 


^res “ Qcon + ^rad + ^vap ^stor ^ucon ^vcon 


( 83 ) 


METHOD OF THE NUMERICAL PROGRAM 


Representations of Physical Properties 

Most of the equations in the previous section contain (at least 
implicitly) quantities representing physical properties of the ablating mate- 
rial or external gas. The density of the ablating material, p, is assumed to 
be constant, but the other pertinent physical properties are considered to be 
temperature dependent. The temperature dependence of the physical properties 
has been left unspecified; however, to obtain numerical solutions, it will be 
specified by formulas with constants that can be arbitrarily chosen and read 
into the computing program. The representations used are: 

Equilibrium vapor pressure of ablating material, atm 

P ve = e w (84) 


Viscosity of ablating material, poises 


B 


q = e 


T-B 


- B 5 


14 


Specific heat of ablating material, cal/g °K 


B, 


c = Be + E X T - ^2 

Thermal conductivity of ablating material, cal/cm sec °K 

Eg o 

K = B g + -f + E 20 T 3 

Average specific heat of external gas, cal/g °K (see eq.. (l4)) 

Cp = E10 + E 3 T w 


( 85 ) 

( 86 ) 

( 87 ) 

( 88 ) 


23 



Transformation of Coordinates and Finite Differencing 

As described in the ANALYSIS AND METHOD OF SOLUTION section, the 
procedure in solving the system of equations presented is to make, in effect, 
all equations auxiliary to equation (l) which is solved with its appropriate 
boundary conditions as given in equations (k6) and (47) or (51)* Equation (l), 
a partial differential equation, is solved by a finite difference scheme; some 
complications are produced by the decrease of length with time between the 
front and back surfaces due to front surface recession. In this work, it was 
elected to keep the number of grid points constant, so the distance between 
grid points was allowed to shrink with decreasing length of the ablator. This 
was accomplished by the following transformation of independent variables from 
t, y to s, t), where 


Then, 


s = t 


_ Y l 
11 L - X(t) 


1- _ 

8y 


L - X(s) 


d_ 

8t) 


a_ 

dt 



( ax/ds ) t) 5 
L - x(s) dr] 


(89) 

(90) 

(9D 

(92) 


This transformation alters somewhat the form of equation (l) and the other 
equations as actually put into the numerical computing program. 


Equation (l) is solved numerically by an explicit (forward difference) 
scheme. In finite difference form, the partial derivatives of the tempera- 
ture, T, are represented as: 



lm,n+i ~ lm,n 
As 


' dTX _ T m+i,n ~ T m-i,n 

driy ~ 2 Ati 

- ‘/m, n 1 



T. 


m+i ,n 


~ 2T m,n + T m-i,n 
(An ) 2 


(93) 

(94) 

(95) 
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Sketch (h) 


where m - nr, m + 1 are grid point 
numbers on the r\ (depth) scale and n , 
n + 1 are numbers on the s (time) 
scale as shown in sketch (b) . Finite 
increments of s and rj are indicated 
by the A symbol. 


Stability and Accuracy of the Finite Difference Equation 

In solving a parabolic partial differential equation by a forward 
difference scheme, there is always a stability requirement to be met. For the 
finite differencing of the transformed version of equation (l), the stability 
requirement turns out to be : 


Z = 


As 

L 



L - X(s)_ 

\ pc / 



(96) 


The stability parameter, Z, is printed out by the computing program for each 
grid point at each time printed. The increments of Ar) remain constant with 
time, but the increments of Ay decrease with time, as indicated by equa- 
tion (90), until ablation is concluded. Thus, Z tends to increase somewhat 
while ablation is proceeding, and this must be considered in selecting initial 
increments. Since Z should not exceed l/2, it is not possible to use the 
computing program to calculate to the point of complete extinction of an 
ablating material. When the present numerical program is used for transparent 
material, the storage limitation requires that the initial value of the ratio 
of length (depth) to the smallest At] be < 1665 (see spacing sketch in 
appendix D ) . 

A gross check on the accuracy of numerical solutions obtained is provided 
by the running energy balance and the cumulative energy balance (see ANALYSIS 
AND METHOD OF SOLUTION section) . An additional check, standard in numerical 
work, can be made by varying As and Ar] and noting the resultant variations 
produced in the numerical solutions. This check, in effect, determines the 
adequacy of representing the derivatives by difference quotients with the 
finite increments as chosen Oeqs . (93) > (9*0 > and (95))- 


Boundary Conditions for Transparent Material 

It was noted in the ANALYSIS AND METHOD OF SOLUTION section that the 
computing program can be arranged to modify the surface energy balance equa- 
tions (46) and (4-7) for the transparent case. Equations (46) and (4-7) are 
rigorous as written, but they require a fine spacing of grid points, and 
therefore small time increments for stability. Particularly near the front 
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face, the derivative of the radiation flux, g, can change by orders of 
magnitude in a short, distance. With a "normal" spacing of grid points, the 
fine structure of the variation of F and g are not well represented, and 
this can result in energy balances that are not accurate (large residuals). 
In the computing program, the front surface energy balance that we actually 
use for the transparent case (Bi S = l) is: 


- V Sy/ = phvVw + + e 3 s(F w - F 1/2 ) ( 97) 

where Fi / 2 is the radiation flux midway between the front surface and the 
first interior finite -differenced grid point (point number 1 + K2 in FORTRAN 
terminology; see spacing sketch in appendix D) . For the back surface energy 
balance. 


( K D bf - %F - ®3 7 (F BF - f BF- 1/ 2 ) ( 96 ) 

where F BF-i/2 is ra diation flux midway between the back surface grid 
point and the nearest interior grid point. The quantities, E 36 and E 37 , are 
constants read into the program and will have values > 0 and < 1. When 
E 36 = E 37 = 0, we recover equations ( 46 ) and (47) exactly. E 36 and E 37 can 
be adjusted to give optimum energy balances; very good energy balances have 
been obtained with E 36 = E 37 = 1. 


ILLUSTRATIVE EXAMPLES 


Examples shown below illustrate the use of the numerical computing 
program as applied to several types of ablation. Calculated and measured 
results are compared for all examples except the last one, which is a Martian 
entry. The disposition of energies for typical examples given is summarized 
in table I. 


Tektite Glass in a Wind Tunnel 

Tektite glasses, ablated at high enthalpies in an arc-jet wind tunnel, 
furnish examples of ablators that both vaporize and melt. Typical comparisons 
between calculated and measured values of surface recession and surface 
brightness temperature are shown in figures 1 and 2. The measured points in 
figure 2 are actually a spread; measurements on other tektite glasses fell 
between these points. The agreement can be seen to be very good, which lends 
confidence that data for flights involving tektite glass can be successfully 
calculated (ref. l) . In both figures the calculations were made for a trans- 
parent glass and for a semitransparent glass, and there is little difference 
between the results of the two methods of computation. The glasses used in 
these examples ablate by melting more than by vaporizing because of the 
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Figure 1. - Comparison of calculated and measured 
surface recession of tektite glass in a wind 
tunnel. 


Figure 2.- Comparison of calculated and measured 
■brightness temperature of tektite glass in a 
wind tunnel. 


moderate enthalpy in the wind tunnel, the low viscosities, and low vapor 
pressures of the glasses. For the glasses in figure 1, about 1 percent of the 
ablation is due to vaporization, and for the glasses in figure 2 , vaporization 
accounts for less than 1 percent. At higher enthalpies the relative amount of 
vaporization increases. Table I gives the disposition of energies calculated 
for the semitransparent glass of figure 1. Because of the small amount of 
vaporization, very little heat is blocked, and most of the incoming energy is 
accounted for by melting. 


Tektite Entry Calculation 


The results of calculating an entry for a typical opaque tektite are 
shown in figures 3 and 4. An entry speed of 11.0 km/sec and an entry angle of 
- 3 O 0 were used for the calculations. These conditions correspond to a typical 
deduced trajectory for a Victoria australite (ref. 2, fig. 22). Figure 3 
shows the calculated values of velocity, surface temperature, surface reces- 
sion, and surface recession due to vaporization. The free -molecule and con- 
tinuum regimes are also distinguished. Time zero is arbitrarily selected as 
"far out" before any appreciable aerodynamic heating has begun. This example 
illustrates the response of a material that vaporizes readily, with about 
14 percent of the ablation due to vaporization and the rest to melting. As 


ABLATION 


BEGINS ENOS 




y. cm 


Figure 3 . - Calculated variation of surface tem- 
perature, velocity, total ablation, and 
vaporized ablation for a Victoria australite 
entering the Earth 1 s atmosphere; R± - 0.8l6 cm. 


Figure U, - Calculated temperature profiles for a 
Victoria australite entering Earth* s atmo- 
sphere; Rj_ = 0 . 8 l 6 cm. 
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the tektite heats up, its surface begins to ablate at a temperature in excess 
of 2000 °K. The surface temperature and ablation rate reach a maximum and 
then fall off as the body slows down. The end of ablation occurs rather 
abruptly, and the remainder of the flight is that of a solid body being aero- 
dynami cally cooled. Measurements of the amount of ablation at the stagnation 
point (refs. 1 and 2) on recovered tektite s yield values not greatly different 
from the 7.3 mm calculated for this example. The calculations indicate that 
for this flight, a negligible portion of the ablation occurred in the free- 
molecule regime. The portion of the ablation in the transitional regime was 
24 percent, compared to the majority of ablation in the continuum regime 
(76 percent) . For smaller tektites and shallow entry angles the percent of 
ablation in the free -molecule and transitional regimes will be greater; for 
large vehicles this portion of ablation is generally small. 

In figure 4 are shown variations of the calculated temperature profiles 
and surface recession at selected values of time. This figure gives a fairly 
complete picture of the internal heating and eventual cooling of the body 
during its flight. The rise in the back temperatures is due to base heating. 
Measurements on recovered bodies by photoelasticity techniques show locked -in 
thermal stresses that vary from a depth of 0.2 to 0.35 cm, corresponding to 
the calculated depth of 0.23 cm (ref. 2). The post-flight observations of 
thermal stresses and deduced surface recession on recovered bodies are com- 
patible with the calculated results. The energy disposition for the flight 
calculated is given in table I. The vaporization that occurs in this flight 
causes substantial heat blockage; the bulk of the energy is accounted for by 
heat blockage, heating and vaporizing, and heating and melt flow, these three 
quantities being of about the same order of magnitude. 


Reentry Flight With Silica Glass Heat Shield 


A calculation was made for a reentry flight of a nose -cone with an opaque 
silica glass heat shield. The vehicle and flight are described by Hidalgo and 

Kadanoff in reference 20. For this trajec- 
tory and this heat -shield material about 
14 percent of the ablation was due to vapori- 
zation which is comparable to the tektite 
entry case previously discussed. The 
recovered reentry vehicle allowed the amount 
of ablation, X, at the stagnation point to be 
determined. The physical -property inputs in 
this case correspond to opaque silica 
(ref. 20), but both the transparent and the 
semitransparent options of the computing 
program were run with the results 



X( transparent) ^ ^ 

X (observed) 


X ( semitransparent ) ^ 

x(observed) ~ ^ 
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The corresponding ratio as calculated 1> y Hidalgo and Kadanoff using their 
quasi-steady ablation analysis was about 1.10 and by Chapman and Larson 
(ref. l) , using an integral method of calculation, 0.92. This illustrates 
that the amount of ablation on simple materials such as glass can be computed 
to the order of 10 - to 15 -percent accuracy by several methods, including two 
of the options of the present computing program. In view of the inevitable 
angle -of -attack variations in flight, which cause the stagnation point of 
maximum heating to wander somewhat over the nose and thus reduce somewhat the 
maximum recession, the observed difference between calculated and measured 
ablation is in the expected direction. The energy disposition for both the 
transparent and semitransparent calculations is shown in table I. It is of 
interest that the two calculations yield energy proportions that are nearly 
the same, although the internal temperature distributions are different 
because of radiative transmission in the transparent case. In both cases the 
total ablation is moderate, so that the actual amount of vaporization is 
moderate and the heat blockage term is relatively small. 


Teflon Model in an Arc -Jet Wind Tunnel 

Under normal ablative conditions, tetraf luoroethylene polymer (Teflon) 
undergoes a surface depolymerization and vaporization of the monomer at a sur- 
face temperature of approximately 760 °K. There is no one specific tempera- 
ture at which the reaction occurs, but a sharply rising reaction rate with 
temperature in this neighborhood essentially controls the surface temperature 
of an ablating model (refs. 21, 22). Under these conditions, the viscosity of 
Teflon remains high, so the process can be said to resemble a sublimation 
(with the reaction rate determined by an Arrhenius type rate equation) . In 
performing the calculations for Teflon ablation, it was assumed that any 
energy involved in possible chemical reactions between the Teflon vapor and 
the external gases could be neglected. 


Comparisons between calculated and experimentally measured surface 
recession for Teflon are shown in figure 5* The experiments by G. Lee and 
R. Sundell (ref. 23 ) were performed in an arc- jet wind tunnel for four values 
of enthalpy. It is seen that the agreement obtained is quite good, with the 
possible exception of the 'JOO Btu/lb total enthalpy case. It is thought that 


l-x*t 



Figure Comparison of calculated and measured 
surface recession of Teflon in wind tunnel. 


the generally satisfactory agreement 
shown in the figure indicates the 
validity of the method calculation and 
also that the physical properties of 
the substance have been adequately 
represented. The front face mass loss 
rate was essentially reaction rate con- 
trolled for the four cases shown in the 
figure . 

The disposition of the calculated 
energies for the 2000 Btu/lb total 
enthalpy case is shown in table I. 

The heat blockage term is fairly large 
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because all the ablated material leaves in the vapor state. The largest term 
is the energy for heating, depolymerizing, and vaporizing. 


Teflon Heat Shield in a Mars Entry 


The calculations for this example illustrate the application of the 
computing program to an entry -with a proposed Mars probe (ref. 24). A spheri- 
cal capsule of 6l.0-cm diameter has been assumed, with a 1-cm thick Teflon 
heat shield, entering the Martian atmosphere in an oriented attitude. Four 
hypothetical atmospheres were assumed as tabulated below. 


Composition (vol.) 

100$ h 2 

91$ n 2 , 9 % co 2 
100$ n 2 

91$ n 2 , 9 $ co 2 


Scale height , km 

7.8 

7.8 

20.0 

20.0 


Subsequent to making these calculations, data obtained from the 1965 Mariner 
occultation experiment have indicated that the scale height of the Martian 
atmosphere is about 9 km* (ref. 25). The assumed atmospheres with the 7*8 km 
scale heights thus appear to be the more realistic ones. Calculations were 
made for an entry velocity of 7*92 km/sec and for two entry angles, - 90 ° and 
-20°. An M/CpA for continuum flow of 3*91 g/cm 2 was assumed for the vehicle, 
and M/A was held constant while the Cp varied through the transition from 
the free -molecule to the continuum regimes. 

The -90° entries have the greater peak heating rates, but the -20° 
entries absorb more total heat and are the more critical from a heat-shield 
standpoint. The heat-shield responses are compared using calculated values 
of the stagnation -point recession as tabulated below. 


Atmosphere 

N 2 small-scale height 
N 2 large-scale height 
N 2 -C0 2 small-scale height 
N 2 -C0 2 large-scale height 
N 2 -C0 2 small-scale height 
N 2 -C0 2 large-scale height 


Entry angle. Total stagnation 

deg point recession, cm 


-90 

-90 

-90 

-90 

-20 

-20 


0.098 

.159 

.121 

.179 

.196 

•341 


Of the -90° entries in the table, the most severe environment would be the 
mixed atmosphere and the large-scale height (although the small-scale height 
appears to be more realistic). For a given scale height, the mixed atmosphere 
gives somewhat more ablation than the nitrogen atmosphere because there is 
more radiation from the carbon dioxide (ref. 17)- The time of exposure to 
heating is roughly proportional to scale height for two otherwise similar tra- 
jectories, and the total heat absorbed is approximately proportional to the 
square root of exposure time (and therefore scale height). This approximate 
relationship between scale height and total recession for a given atmospheric 
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composition can be deduced from the table. The strong dependence of the 
heat-shield response on any uncertainty of knowledge of atmospheric scale 
height is of interest, since this trend will presumably apply to any heat- 
shield material and any planet. 

The disposition of energy for the -20° entry with the mixed atmosphere 
and short -scale height is shown in table I. As with the wind-tunnel results 
for Teflon, the two large energy terms are the heat blockage term and the 
term that accounts for heating, depolymerizing, and vaporizing. In the envi- 
ronment of the -20° Martian entry, the considerable amount of ablation of 
material to the vapor state accounts for the very large heat blockage term. 


MARS ENTRY ASSUMED ATMOSPHERE 

ENTRY VELOCITY = 7.92 km/sec 91% N 2 , 9% C0 2 (VOL) 

ENTRY ANGLE = -20* SCALE HEIGHT = 20.0 km 

TOTAL MATERIAL LOSS = 1.56 kg = 3.43 lb 



Figure 6. - Approximate calculation of response 
of a Teflon heat shield in a Mars entry. 


As an illustration, the analysis 
has also been used in an approximate 
manner to calculate the quantities of 
interest around the front hemisphere of 
the spherical capsule for the -20° 
entry in the mixed atmosphere with the 
(probably overly severe) large scale 
height. The results are summarized in 
figure 6 which shows the variation of 
total recession and front and back face 
maximum temperatures around the hemi- 
spheric heat shield as well as the 
total mass loss for this hypothetical 
case. These results illustrate how 
variable material thickness may be used 
in heat-shield design. 


TABLE I.- TYPICAL ENERGY BALANCES (percentages) 



Tektite glass 
( semitransparent) 
wind tunnel 
(fig- 1) 

Tektite 
Earth entry 
(opaque) 
(figs. 3,4) 

Silica glass 
heat -shield 
Earth entry 
(transparent) 

Silica glass 
heat -shield 
Earth entry 
(semitransparent) 

Teflon wind 
tunnel (fig. 5) 
h s =2000 Btu/lb 
p s =0.34 atm 

Teflon 
Mars entry 

n=-2o° 

v ooi=7-92 km/sec 

Convection in 
(hot wall) 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

Heat blocked 

2.9 

31.9 

11.4 

11.5 

35-4 

71.9 

(Wet convection in) 

(97.1) 

(68.1) 

(88.6) 

(88.5) 

(64.6) 

( 28 . 1 ) 

Net radiation in 

- 13.3 

-TO .8 

-28.4 

- 27.6 

-4.8 

2.4 

Heating and 
vaporizing 

5.2 

24.9 

21.5 

21.2 

50.1 

25-5 

Heating and melting 

69.1 

27.2 

28.6 

31-6 

0 

0 

Stored 

4.6 

6.3 

10.4 

6.9 

9-9 

5.1 

Error 

4.9 

-1.1 

-0-3 

1.2 

-0.2 

-0.1 
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CONCLUDING REMARKS 


A generalized analysis of stagnation-point ablation has been presented 
for solving a variety of problems involving melting and vaporizing, subliming, 
or surface chemical reactions. The flexibility of the analysis has been 
demonstrated through the presentation of several varied illustrative examples. 
In general, it is expected that accuracy of answers obtained will depend 
largely on the degree of knowledge of the physical, chemical, and thermody- 
namic properties of the ablating material, as these are necessary inputs for 
the computing program. The procedure of relating calculations for a given 
material to experiment wherever possible lends confidence to calculations for 
the same material exposed to other conditions which cannot be verified by 
observation. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., July 22, 19 66 
129-03-12-01-00 
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APPENDIX A 


PRINCIPAL NOMENCLATURE 


In performing computing machine calculations, some purely FORTRAN 
quantities are used, particularly among input data, which have no counterpart 
among the symbols listed below. These quantities are in appendix D, wherein 
all FORTRAN quantities are listed. 

a y direction body force per unit mass, cm/sec 2 (acceleration in 

flight) 


*Dg 




A C qjA cv , 

A cm 


A 1 

a cv 


Ax 

A 2 

A3 

A 4 

A 2D 

b 

B 

B 2 

b 3 

b 4 

b 5 


component of deceleration due to aerodynamic drag normalized by 
gravitational acceleration of planet , dimensionless (eq. (69)) 

absolute value of aerodynamic acceleration normalized by gravita- 
tional acceleration of planets dimensionless (eq. (70)) 

frontal area, cm 2 

free -molecule accommodation coefficients for heat transfer, mass 
loss, x momentum (for surface shear), respectively, dimension- 
less 

corrected mass -loss accommodation coefficient (eq. (39c)), 
dimensionless 

constant, defined by equation (27) 
constant, defined by equation (ll) 
constant, defined by equation (32a) 
constant, defined by equation (15) 

melt-off parameter, defined by equation (65) > dimensionless 
Sutherland constant, °K (eq. (C38)) 

Arrhenius frequency factor, cm/sec (eq. ( 39 e)) 
constant in vapor pressure (eq. ( 84 -)) 
constant in vapor pressure (eq. (8^)) 
constant in viscosity (eq. (85)) 
constant in viscosity (eq. (85)) 
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By 

b 9 

Bn 

B14 
B 16 

c 

c 

c p 

c p 

Cd^°dcA 

C DFM / 

Ci 
C 2 
C 3 
c 4 
C 6 

D 

e aZ 

e rr 

Eiji 

E 2 (z) 

Ei 

E2 


constant in specific heat (eq. (86)) 
constant in specific heat (eq. (86) ) 
constant in thermal conductivity (eq. (87)) 

constant in convective heat blockage factor (eqs. (28), (30)) 

(See discussion following eq. (C 7 )*) 

constant in viscosity (eq. (85)) 

constant; B 1S = 1.0 for transparent case; B 16 = 0 for opaque and 
semitransparent cases (eqs. ( 46 ), ( 47 ) j (51)) 

specific heat of body material, cal/g °K 

heat capacity per unit area of backing material, cal/cm 2 °K 
(eq. ( 51 )) 

specific heat of a gas at constant pressure, cal/g °K 

average specific heat, external gas, cal/g °K; defined by equation (l 4 ) 

drag coefficient, continuum drag coefficient, free -molecule drag 
coefficient, respectively, dimensionless 

constant in nose radius (eq. (60)) 

constant in nose radius (eq. (60)) 

constant in M/A (eq. (56)) 

constant in M/A (eq. (56)) 

constant, vorticity correction in equation (15) 
free -stream density, g/m 3 

allowable error in T w (selected), °K; allowable disagreement 
between T w obtained from equations (l) and ( 46 ) 

error in T v after last iteration, °K; will be < e a ^ 

Arrhenius activation temperature, °K (eq. ( 39 ^)) 

exponential integral (second degree), defined in equation ( 43 ) 

constant in specific heat (eq. (86)) 

constant in thermal conductivity (eq. (87)) 
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e 3 

e 4 

e 5 

e 6 

e 7 

E s 
Eg 
Eio 
En 
E 12 

Ei3 

E14 

El 5 

El6 

E17 

E 18 

E19 

E 2 o 

E33 

E35 

E 3 6^E37 

E3S 


constant in average specific heat (eq. (88)) 
constant in gas -cap radiation (eq. (44)) 
constant in gas-cap radiation (eq. ( 44 )) 
constant in gas-cap radiation (eq. (44)) 

constant to account for shift of vaporization equilibrium 
(eq.. (26)) 

constant in expression for shear blowing factor (eq. (34)) 

constant, defined by equation (58) 

constant in average specific heat (eq. (88)) 

constant in expression for p 21 (eq. (6l)) 

constant in expression for p 21 (eq. (6l)) 

constant, defined by equation (50) 

constant depending on body shape and flight conditions in drag 
bridging (eq. ( 57 )); see also equation (C 46 ) 

constant accounting for body force in wind tunnel expression for 
P ,f (eq. (lOa) ) 

constant used in tumbling correction (eq. (12)) 
constant in nose radius (eq. (60)) 
constant in M/A (eq. (56)) 

constant, defined by equation ( 49 ) as the ratio of turbulent to 
laminar base heating 

constant in thermal conductivity (eq. (87)) 

constant used in expression for front face emissivity for semi- 
transparent body (eq. ( 45 )) 

constant, asymptotic value of \|r (eq. (28)) 

(See discussion following eq. (C 7 ).) 

constants used in equations (97) and (98) to modify surface energy 
balances for the transparent case 

average ambient enthalpy for flight case, km 2 /sec 2 (eq. (55)) 
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F 

F RS 

F l/2 


F BF-l/2 


F 

F t 

g 

S P 

h 

hs 

hv 

hoo 

h 

h 

K 

K A 

Ka 

K 


Ktu 


L 

Le 
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ratio of actual viscosity to undissociated (Sutherland) value, 
dimensionless 

radiation flux, cal/cm 2 sec (eq. (ij-l)) 

reradiation rate from front surface, cal/cm 2 sec (eqs. (71), (72)) 

radiation flux midway between front surface and first interior 
finite -differenced grig point, point number 1 + k2 in FORTRAN 
terminology; see Spacing Sketch in appendix D; cal/cm 2 sec 

(eq. (97)) 

radiation flux midway between back. surface grid point and nearest 
interior grid point, cal/cm 2 sec (eq. (98)) 

ratio of surface recession rate due to vaporization to total sur- 
face recession rate, dimensionless (eq. (66)) 

ratio of surface recession due to vaporization to total surface 
recession, dimensionless (eq. (67)) 

gradient of radiative flux, dF/dy, cal/cm 3 sec 

gravitational acceleration of planet, cm/sec 2 (eq. (53^0) 

enthalpy, cal/g 

stagnation enthalpy, cal/g 

latent heat of vaporization, cal/g 

average ambient enthalpy for flight case, cal/g (eq. (55)) 

h/h s , dimensionless 

(l - h)/(l - 3 %), dimensionless 

thermal conductivity, cal/cm sec °K 

mean free path constant, moles (eq. (C 40 )) 

Reynolds analogy factor (eq. (31))^ dimensionless 

fractional time center point is exposed to stagnation conditions, 
dimensionless (eq. (l2)) 

correction factor for heat transfer and other quantities due to 
oscillation in flight, dimensionless (eq. (12)) 

initial depth of material, cm 

Lewis number, dimensionless 



L/D r 

m 

• 

m 

M 

M oo 

M 

n 

P 

Pt 2 

Pv 

*Ve 

Pvm 

P 

p" 

q. 

^oc 

q FM 

% 

9oc 

vort=o 

% 


lift/drag ratio , dimensionless 

molecular weight 

mass loss rate, g/cm 2 sec 

mass of body, g; also grid point number in finite difference 
computation 

free-stream Mach number,, dimensionless 
me/m v , dimensionless (eq. (29) ) 

index of refraction, dimensionless; also exponent in equation (30) 

pressure, dynes/cm 2 or atm, as specified 

pressure downstream of normal shock, atm 

actual vapor pressure, atm 

equilibrium vapor pressure, atm 

modified equilibrium vapor pressure, atm (eq. (26)) 

ratio of pressure of external gas to modified equilibrium vapor 
pressure of ablated vapor, dimensionless (eq. (26)) 

negative of second derivative with respect to x of external 

pressure with a correction for body force, dynes/cm 4 ; defined in 
equation (8); evaluated in equation (10) 

heat-transfer rate, cal/cm 2 sec; also dynamic pressure, dynes/cm 2 

surface convective (continuum) heat -transfer rate with no blowing, 
cal/cm 2 sec (eq. (15)) 

surface convective (free molecule) heat -transfer rate, cal/cm 2 sec 
(eq. (17)) 

surface convective heat -transfer rate with no blowing, bridged 
between q oc and q cal/cm 2 sec (eq. (20)) 

q OQ corrected for tumbling or oscillation, cal/cm 2 sec (eq. (21)) 

q oc without vorticity, cal/cm 2 sec; C 6 = 0 in equation (15) 

gas-cap radiation rate, cal/cm 2 sec (eq. (44)) 


# 
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surface convective (continuum) heat- transfer rate with blowing , 
cal/cm 2 sec (eq. (l6)) 

CLlJrW 

surface convective heat -transfer rate (all regimes), bridged 
between and q^, cal/cm 2 sec (eq. (l8) or (22b)) 


q^ v corrected for tumbling or oscillation, cal/cm 2 sec 
(eq, (19) or (23)) 

9 -wRi 

initial conib Ined convective and radiative heating rate, 
cal/cm 2 sec (eq. (B3a)) 

Energy transfer rates listed below pertain to an energy balance, and, 
where applicable, are combined rates for front and back surfaces. 

q 

u con 

total convective heat- transfer rate into material, cal/cm 2 sec 
(eq. (75)) 

^rad 

net radiative heating rate into material, cal/cm 2 sec (eq. (76)) 

^vap 

energy rate due to vaporization (positive if energy is released 
into the material) , cal/cm 2 sec (eq. (77) ) 

^stor 

rate of increase of stored energy In the material, cal/cm 2 sec 
(eq. ( 78 )) 

^•ucon 

x direction convection energy rate of the material (positive 
out), cal/cm 2 sec (eq. (79)) 

q 

^-vcon 

y direction convection energy rate of the material (positive 
out), cal/cm 2 sec (eq. (80) ) 

^•res 

residual in energy rate balance, cal/cm 2 sec (eq. (8l) ) 

^con'^rad, 

Qvap^^stor^ 

^ucon^^vcon 

J 

time integrals of corresponding q values, cal/cm 2 (eq. (82)); 
terms in total energy balance 

Qres 

residual in total energy balance, cal/cm 2 (eq. (83)) 

R 

nose radius, cm 

R b 

body radius, cm (eqs. ( 48 ), (50)) 

%/F 

ratio of base to front -face laminar convective heating 

normalized by respective enthalpy potentials, for exposed 
back surface, dimensionless (eq. ( 48 )) 

Re 

Reynolds number 
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~R e ff effective coefficient of reflection for planar radiation, 

dimensionless (eq. (42)) 

Re s Reynolds number based on enthalpy velocity (eq. (C20) ) 

Rg universal gas constant, ergs/mole °K 

Rp planet radius, km (eq. (53)) 

s transformed time coordinate, sec (eq. ( 89 )) 

S collision cross-section area, cm 2 (eq. (C40) ) 

atmospheric scale height, km (eq. (54)) 

S^. initial scale height for entry into arbitrary atmosphere, km 

(eq. (B6) ) , and contained in equation (BIO) 

Sh 1 ^Sh 2 ^Sh 3 successive scale heights in arbitrary atmosphere, km: 


s h = s hi 

when 

Poo — Poo 2 

s h “ s h 2 

when 

Po0 2 ^ Poo 

Sh = s h3 

when 

Poo > P003 


time, sec 

initial time, sec; time at which front face temperature, T w , 
arrives at assumed for wind tunnel cases (eq. (B4)) 

temperature, °K 

brightness temperature (emissivity unity) , °K 
reference temperature, °K 

viscosity thickness temperature, °K; temperature at which 
M- = e.fXy. (at depth A^) 

horizontal component of trajectory velocity, km/sec 

x direction velocity of external gas at edge of boundary layer, 
km/sec (eq. (ll) ) 

velocity of material in x direction, cm/sec 
u .10 5 , cm/sec 


9 sec 
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u* 


V 

V 


v ( 

V 

V 


sr 


x 


X r 

y 

Y 


Z 


a 

as 

7 

& 

8X 

8 * 

A 






€ 


1 

?\ 


u/u e , dimensionless (eg.. (C2l) ) 

vertical component of trajectory velocity (positive upward), km/sec 
velocity of material in y direction, cm/sec 
surface recession rate, cm/sec (eq. (6a)) 

enthalpy velocity, km/sec; defined as V 2 = O.OO836 hs (eq. (13)) 

free -stream (trajectory) velocity, km/sec (eq. (52a)) 

longitudinal coordinate along meridian, cm 

flow line ratio, dimensionless; defined by equation (68) 

transverse coordinate normal to surface (inward), cm 

boundary-layer transverse coordinate, cm; appendix C 

stability parameter for finite differencing, dimensionless; 

(eq. (96)); must be < l/2 

absorption coefficient, internal radiation, cm -1 (eq. (tl)) 
absorption coefficient, gas-cap radiation, cm -1 (eq. (tl)) 
trajectory angle, deg, positive above horizontal (eq. (52b)) 
increment (appendix C) 

mass fraction, dimensionless (eq. (CIO)) 
displacement thickness, cm 

thermal thickness, cm (eq. (63)); also increment e.g.. 

Ah = enthalpy potential, cal/g 

unspecified characteristic boundary- layer thicknesses, cm 
(appendix C) 

viscosity thickness, cm; depth at which q = e.q w 

surface emissivity, opaque or semitransparent; effective emissivity, 
transparent (eqs. (73) > (7*0 )> dimensionless 

transformed y coordinate, cm (eq. (90 )); also dummy variable 

mean free path, cm 


M- 
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viscosity, poise 



p density of ablating material (constant), g/cm 3 

p density ratio across normal shock, dimensionless; for -wind tunnel 

" 1 cases, assigned; for flight cases, equations ( 6 l) and ( 62 ) 

p^ free -stream atmospheric density in Earth sea-level atmospheres, 

dimensionless, D/1226 

Poo >Poo values assigned to p^ at which changes of scale height in an 
2 3 arbitrary atmosphere occur; see S^, S]^ 

a Stefan constant, 1.369X10" 12 cal/cm 2 sec °K 4 ; also Prandtl number, 

dimensionless 

t shear, dynes/cm 2 

cp approximation of the ratio of stored energy to the stored energy 

associated with an exponential temperature profile, dimensionless 

(eq.. (64)) 

X surface recession, cm 

X va p surface recession due to vaporization, cm; X vap = F^X (eq. ( 67 )) 

Xi characteristic recession depth, cm, used in tumbling correction 

(eq. ( 12 )) 

\|/ convective heat blockage factor, dimensionless (eqs. (l 6 ), ( 28 )) 

\| r modified convective heat blockage factor, dimensionless (eqs. ( 22 ), 

(24)) 


a 

BF 

c 

cha 

d 

e 

eq 

FF 


Subscripts 

actual 
back face 
continuum 

change of wind tunnel conditions 
diffusion 

external gas or outer edge of boundary layer 

equilibrium 

front face 
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FM free molecule 

i initial; this subscript can be combined with the others 

max maximum 

0 no blowing 

oc no blowing, continuum 

off shut off of wind tunnel 

r reverse 

ref reference 

s stagnation (or settling chamber) 

u undissociated 

v vapor expelled 

w wall (front face) 

wc wall, continuum 

wd wall, diffusion 

wFM wall, free molecule 

1 dummy variable 

2 behind normal shock; also average condition between shock wave and 

body (appendix C) 

00 free stream 

Superscript 

* x derivative 




APPENDIX B 


EQUATIONS FOR STARTING- VALUES 

To start the solution to equation (l), it is necessary to assign initial 
conditions. These will normally consist of a relatively low temperature pro- 
file which can exist before the onset of ablation. The particular selection 
of initial conditions is generally not critical as their influence damps out 
in a short time. The initial temperature profile that we assume is an expo- 
nential type . 

Ti(y) - T 0 + (T wi - T 0 ) 


-y/A ± 




“^BF 


(Bl) 


If T wi is selected near T 0 , this profile amounts to a small perturbation 
on the constant T 0 profile. We can take the y derivative of Tj_ at 
y = 0, equate it to the ratio of the heat flux (eq_. (46)) to thermal conduc- 
tivity at the wall* and solve for the initial thermal thickness, A^. 


•^wi^wi ^o^ 
'iwi + W “ B i6)q.Ri 



(B2a) 


We use the more simple approximation: 


A ^vi^wi l^o) 

1 ~ %ri + Cl - B ls )q Ri 


(B2b) 


where Bi© = 1 for the transparent case and B 16 = 0 for the opaque and semi- 
transparent cases. As described below in the section. Flight Cases, one 
option allows A± to be assigned a value instead of obtaining it from equa- 
tion (B2b). The initial convective and radiative heat fluxes needed in equa- 
tion (B2b) are obtained differently for wind tunnel or flight as described 
below. In calculating q vi and q Ri , the initial free-stream density, D^, is 
needed. For wind tunnel calculations, Di will be known; for flight, Dj_ can 
be assigned or calculated as described below. 


Ablation in a Wind Tunnel 

The starting of a wind tunnel is visualized as a sudden step of a heat 
flux. We define a combined initial heat flux, q w Rj_, as the sum of the initial 
convective and radiative fluxes. 


q wRi " q wi + q Ri 


(B3a) 
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We assume that ^ = 1 and that the convective heat transfer is in the 
continuum regime (OpMi no ^ calculated). With q^p = q^^ we have: 

* 

<lwRi = 'loci + ^Ri (B3b) 

We evaluate ( l oc ± from equation (15) and q^-p from equation (44), with 
Kj_ u = 1, using T-^-i for T-^-. If T w i is small, q Q ci an( ^ therefore 
are approximately constant for a short period of time (eq. ( 15) ) - The classi- 
cal conduction problem with a constant heat flux and constant properties 
(ref. 2 6, p. 56) can be used as an approximation for this case to determine 
the time at which the front face temperature arrives at the assumed Twi- 
This turns out to be: 


t i = 


( T wj T Q ) rt P SrL^wi 

^wRi 2 


(B4) 


One can set T v p > T Q ; the greater value is not necessary, but it gives the 
computing program a smooth start. 


Flight With Arbitrary Initial Conditions 

In starting the flight calculations, we will use an assigned initial 
velocity and a flight path angle, 7^, at (arbitrary) time, t^ = 0. The 

initial atmospheric density, Dj_ (equivalent to a starting altitude), can be 
assigned, as well as an initial thermal thickness, Aj_. With an assumed Twi; 
the initial profile is determined from equation (Bl) . Using this starting 
procedure, we do not require that Aj_ be consistent with the relationship 
given in equation (B2b) . However, we calculate q w j_, q Ri , and q v ^p as fig- 
ures of interest as described below in the section. Initial Conditions for 
Entry Flight . 


Initial Conditions for Entry Flight 

An alternative starting procedure, valid for an entry ’flight, is to use 
an assigned V ^ and 7^ and an assumed. Twi; and to calculate the entry into 

an exponential atmosphere [D = c e -(Alt/Sh-p) ] which will raise the front face 
temperature to the assumed T^i- The convective heating during this initial 
part of the entry will be considered to be of the free -molecule type (and we 
do not calculate q^p). can write (eq. (19)): 


^-fwi^u 


where 

^Jfwi ~ ^FMi 
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The evaluation of K^ u in equation ( 12 ) specializes to 

K t u - K 


so we have 


q wi “ Kq PMi 

We evaluate q^,^ from equation (17), and using equations (B 3 a) and ( 44 ), we 
have 


KAcq^iVi 


Ebv E s 


q wRi = o.683g- < V i 2 - °- 008 3 6 = P Twi) + KE^iDi^i 


(B5) 


where we use Vj_ as an approximation for V 


The quantity l^-Ri is heating rate at time t^ = 0 . Up to time 

hi = 0, the heat flux will be approximately an exponential function of time, 
through the exponential variation of D. This is similar to a classical prob- 
lem in heat conduction (ref. 2 6, p. 45) which yields an exponential tempera- 
ture profile. In equation (B 5 ) we can replace B± by D as an exponential 
function of altitude; we can approximate the slightly varying enthalpy poten- 
tial as a constant (the value in eq. (B 5 )); and we can integrate the heating 
flux over time from t = -°q to t = t^ = 0 . We obtain, then, the total heat 
absorbed: 


^total 


Sh-^KA. C qD (V-^ - O.OO836 CpT-^q) 

+ 

O.O836 | sin 7-jJ 


— Rg Eq 1 

ShjKBjgigi Vj 

Eg | sin 7 ± j 


(B6) 


We can also approximate the total heat absorbed as 


^total 



T 0 )dy 


(B7) 


When we substitute the profile equation (Bl) into equation (B 7 ), we obtain 

^total ** P c wi ^wi - ^o^i^ + %5 /F ^ (^ e ^ (B8a) 

or, approximated further, 

^total ^ pCw i ^wi “ ^o^i (B8b) 

We substitute q wi and q Ri (the first and second terms, respectively, on the 
right side of eq. (B 5 ) ) into equation (B 2 b) for A-^, put this into equa- 
tion (B8b), and have finally 
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(B9) 


We now eliminate 
tion of the form 


Qtotai ^ rom equations (b 6) and (B9) and we obtain an equa- 





+ 



(B10) 


where Ki and K 2 are constants and D-j_ is the only unknown. 


The procedure, then, for starting the entry calculations is as follows. 

We assume a T vi. and we know the scale height, Shj_, far out in an atmosphere. 
We use equation (BIO) to calculate we can then calculate q wi and q^j_ as 

the first and second terms, respectively, on the right side of equation (B5) • 
Finally, we obtain from equation (B2b), and we put A^ into the profile 

equation (Bl) . We can see that the assumption of ^wi fixes the D-j_, or, in 
effect, fixes the altitude at which we start time zero. For this case, 

T^i > T 0 is necessary in order to have a finite starting altitude. 
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APPENDIX C 


BRIDGING BETWEEN FREE -MOLECULE AND CONTINUUM REGIMES 


For a number of applications or situations, it is known that the 
transitional regime between free -molecule and continuum flow must be consid- 
ered. The bridging formulas used in this work for the transitional regime are 
presented without derivations, in* the ANALYSIS AND METHOD OF SOLUTION section 
as equations (l 8 ), (20), (24), ( 37 ), (40), and (57)- The derivations of these 
equations, based principally on a simple kinetic theory model, are given in 
this appendix. These bridging formulas replace previously used equations 
which were of purely empirical form (ref. l) . 


Front Face Normal Velocity 

For the front face normal velocity or mass loss rate, a bridging 
relationship is required between the free-molecule and continuum regimes. We 
will consider first the free-molecule or reaction rate -controlled regime, and 
we will distinguish between the evaporation or sublimation case and the 
chemical reaction case . 

For the case involving evaporation or sublimation, we write the Langmuir 
equation (eq. ( 76 ) of ref. 27 ) for the mass loss rate into a vacuum as: 


/ T m v 

®FM = A cvP ve C dy J 2rt R g T w ^ Cl) 

where the constant, C^y, is the pressure of a standard atmosphere in 
dynes/cm 2 so that p ve is measured in atmospheres, m v is the molecular 
weight in g, and Rg is the universal gas constant in erg/mole °K. Equa- 
tion (Cl) is based on the rate of molecules crossing a unit area, in this 
case impacting against a unit surface area, and the accommodation coefficient, 
A cv , is the fraction of the molecules that stick to the surface and condense. 
At equilibrium the evaporation and condensation rates are equal so there is 
then no net evaporation rate. At a given temperature, the rate of surface 
impacts, and therefore the rate of condensation, is proportional to the actual 
vapor pressure above the surface, p . So we can write for the condensation 
rate : 

p 

m r = ( C2 ) 

v ve 

and for the net rate of evaporation (or sublimation) 

m = - m r (C3a) 
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If the equilibrium vapor pressure is shifted or modified by the presence of 
other gaseous materials in the boundary layer, then the largest actual vapor 
pressure at the surface, p v , that can be reached is the modified equilibrium 
vapor pressure, p w (see eq. (2 6)), rather than p ve . We accordingly modify 
equation (C 3 b) and use the approximate expression: 


m = m FM ( 1 - 



For the chemical reaction case we rewrite the Arrhenius rate equation 
for a first order reaction (eq. (39e)) as 


m FM = P^ e 


-Et/T, 


w 


(C 4 ) 


When the modified equilibrium vapor pressure, p^, exists above the surface, 
the net reaction rate is zero, or the reverse reaction rate equals the forward 
reaction rate. At a given temperature, the rate of impact of vapor molecules 
with the surface will again be proportional to the actual vapor pressure at 
the surface, p v . We can assume that the reverse reaction rate, is propor- 

tional to the rate of surface impacts and therefore to the actual vapor pres- 
sure, p v , and we again have equation (C 3 c) for the net reaction rate. 

We now consider the continuum or diffusion controlled regime for the 
front face mass loss rate. The limiting value of the diffusion controlled 
front face normal velocity, |v wc |, is given in equation (38), and we define 
% = p|v vc | (obtained by putting p on the left side of eq. (38)), as the 
limiting value or maximum mass diffusion rate. As noted in the discussion 
following equation (38), this evaluation is obtained using P and \|/ which 
ultimately depend on the modified equilibrium vapor pressure at the surface, 
p^. As distinct from the theoretical or limiting value, we next consider 
the actual mass transfer rate by diffusion. The actual rate can be shown to 
be approximately proportional to the actual vapor pressure at the surface, 

P v , by inserting actual, rather than limiting, evaluations of P and into 
equation (38) . In these actual evaluations, P will be based on p v rather 
than p m and the \|r asymptote is assumed to be ^ 0 with B n w 1 . We then 
obtain 



This is straightforward when p m < p t . When p m > Pt 2 j a modification is 
necessary because in equation (38) we have, in effect, given p m an upper 
limit of 0 . 999 j 999 Pt 2 * case we have 


Let 


m 






(C 5 b) 


“d 


t 



and we have equation (C5a) with m^ 1 replacing m^. 


(c6) 


If a quasi-steady state is assumed, there is no accumulation or depletion 
of vapor in the boundary layer. Thus the actual mass loss rate calculated by 
the free -molecule (reaction-rate) method, and the actual mass loss rate calcu- 
lated as a diffusion process must be the same; physically this means that 
the actual vapor pressure at the surface, p v , must have a value such that ST 
obtained from equations (C3c) and (C 5 ) will be the same. We eliminate Pv/Pvm 
between equations (C 3 c) and (C 5 a) and obtain 


1 


m 



(C 7 ) 


At high mass loss rates, one can surmise that diffusion control may merge into 
a hydrodynamic control with an interface between the two fluids. In this case, 
equation (38) may yield values of the ’’diffusion” rate that are too large at 
the high mass loss rates. Finite values of the ”dif fusion” rate will be 
obtained from equation (38) when the \|r asymptote E 35 = 0 , and Bn = 1 . 

For the (somewhat unusual) situation when P-vm/Pt^ >0.8, it is suggested 
that one use E 35 = 0 and Bn & 1 ; this insures that = p|v wc |, as calcu- 
lated using equation (38), does not become unrealistically large. ^ For the 
condition with Pyrri > Pt p j we should actually replace m^ with mq* in 
equation (C7). However, this situation will generally be one of rate control 
in which will be small relative to r% (whose calculated value may be 

too large) and % f will be still larger. As an approximation for all 
conditions, we will accordingly, use equation (C7) with M3. 

When we cancel out the constant density, p, from equation (C 7 )> we have 
equation ( 40 ) for the front face normal velocity, v w . A somewhat similar line 
of reasoning for the evaporation case is in reference 28 although bridging 
equations are not presented. These bridging forms (eqs. ( 40 ) and (C7)) are 
considered to be valid approximations over the complete spectrum from rate 
control to diffusion control; the use of the bridging relationship automati- 
cally places control in the proper regime. Some verification of this bridging 
relationship is given in figure 5 by the comparisons between calculations and 
experiment for the ablation of Teflon in a wind tunnel. These examples are 
essentially reaction rate controlled, except for h s = 3000 Btu/lb which is 
considered to be in the transitional regime. 
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Convective Heat Transfer 


To calculate heat transfer in the transitional regime, we require a 
bridging between the continuum and free -molecule heat transfer values given by 
equations (l 6) and (17). The bridged result is shown in equation (l8). A 
very simple first-collision model is used in the analysis. A typical packet 
of free molecules is assumed to enter the boundary layer and make a first 
collision with molecules already there. The free molecules then become part 
of the continuum boundary layer with average energy equal to that of the 
boundary layer at the point of the collision. The energy given up by the free 
molecules on collision is assumed to be ultimately taken up by the wall (by 
successive collisions in the continuum boundary layer and impact with the 
wall), since there is no piling up of energy in a quasi-steady state boundary 
layer. (Some of the free molecules will make their first collision with the 
wall and give up energy directly. ) 


We use a normalized enthalpy, h = 
in sketch (d). We define an effective 



entering the boundary layer and making 
allows us to write: 


h/h s , and the coordinate system shown 
collision thickness, Aq, such that 
collisions occurring within Aq have 
an effect on wall heat transfer, while 
collisions occurring outside of Aq 
have a negligible effect on wall heat 
transfer; Aq is thus a kind of bound- 
ary layer thickness. We can rewrite 
equation (17) as 

Q.FM = Ko/DVjhsU - h w ) (C8) 

where Ki contains the conversion of 
units and the accommodation coefficient, 
A C q (assumed to be approximately unity). 
Then, for a small packet of molecules 
its first collision, our assumed model 


&q.FM = [h s Ki(DVj]SX(l - h) (C9) 

where 5X is the mass fraction of molecules that make the first collision 
between Y and Y + dY (or between h and h + dh) in the boundary layer; bqp^ 
is heat given up by the fraction SX. In using K x in equation (C9), we are 
assuming that the fraction of energy given up in the molecular collisions is 
the same as the fraction given up by collisions with the wall (approximately 
unity) ; this is thought to be within the framework of approximations being 
made in the derivation. The mass fraction, 5X, is evaluated in terms of the 
mean free path, A (ref. 29, eq. (103-7) )• 

6X = e" Y A f (CIO) 
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The mean free path. A, is actually a function of Y since A varies inversely 
with density in the boundary layer. We will consider A to have some con- 
stant mean value in the boundary layer in order to perform an integration 
(below) . This assumption appears to be within the framework of approximations 
being made. In accord with our quasi- steady-state assumption, we write 


< fyw = 5q -FM 

all 5X 


and using equation (C9) we have 


Q'ty 


w 


all 5X 


[hsK^DVj ] [hsKxCDVj] 

¥e insert SX from equation (CIO) and have: 

Srw 


5q -FM 

= ) 5X(l - h) 

ur M l — i 


all 5X 


[h s K ± (mJ] 


(1 - h)e 


x 


all SX 


We will use : 


= - 3^4- - f (£ 

1 - h w W 


and we can write: 




w 


[hgK^DVj] 


- (l - hy) 


fi * « V 

A 


(C11) 


(C12) 


(C13) 


(C14) 


(C15) 


The last term in the "bracket accounts for heat transfer from the molecules 
whose free path is greater than Aq (ref. 29, eq. ( 103 - 8 )). Using equa- 
tions (C 8 and C14) we have finally 


V 


w 




•FM 


^ „- Y A e + »-V A 


A 


(Cl 6 ) 


We will use the simplest form for f(Y/Aq) in equation (Cl4) (see sketch (d)). 

h = fiVf (C17) 

^q 
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We perform the integration in equation (Cl 6 ) and obtain 


V 


= In 


-PM 


7±_ 

K 



(CIS) 


To consider a variation in the ratio, A/Aq, we can visualize a change in 
ambient density, D. The mean free path, A, varies as D~ x , and Aq (being a 
kind of boundary layer thickness) varies as D~ 1/2 , so the ratio, A/Aq, varies 
as D” 1/2 . In the free-molecule limit, A/Aq becomes very large, and in the 
continuum limit, very small. When A/Aq is large, q^ w approaches in 

equation (Cl 8 ), as it should. When A/Aq is small, q^ w must approach q^ c . 
Using equation (Cl 8 ) -we can write for small A/Aq 

<tyv = V = ^ 1FM (C19a) 


_ qpM 
^ = V 

The quantities, q^ and q^ c , are calculated values for given (the same) condi- 
tions. Since the characteristic thickness, Aq, has not been specified, it can 
now be given the value that satisfies equation (C19b) for the conditions 
imposed. The mean free path. A, has been considered to be some mean value in 
the boundary layer, and it can be assigned a convenient value, say A 2 , with 
Aq required to satisfy equation (C19b) with this A. It can be shown that 
the form of equation (C19b) is a consistent relationship by making use of 
equations (5-5), (6.23), (6.25), ( 6 . 56 ) of reference 30, and relating A to 
a characteristic thickness (displacement thickness, 5*) . 

We substitute equation (C19b) into (Cl 8 ) and we have equation (l 8 ) as 
the bridging equation for convection heat transfer. 



V * V ( q - e' q M /q ^) ( 18 ) 

Equation (l8) has the form of a monotonic function of q^/q^ c , and it has 

the correct free-molecule and continuum regime asymptotes. The derivation 
has been essentially performed from the free-molecule end, and the agreement 
with the continuum end has been forced. In the derivation, functional forms 
for f(y/Aq) other than the simple one chosen (eq. (CI 7 )) can be used. The 
result will be similarly behaved, but more complicated, bridging equations. 

It is possible that the more complicated equations that can be obtained may 
give better agreement with measured data for some specific types of heat 
transfer . 

The alternate forms of the bridging relations (eqs. (20) to (24)) are 
obtained directly from equation (l8) as outlined in the ANALYSIS AND METHOD OF 
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SOLUTION section. The alternate forms yield the quantities , q^, q , and \|f 
(which are mainly of conceptual interest) . " 


Comparisons of calculated bridging with experiment are available for the 
case with no blowing. Experimental data, obtained from figure 5 of refer- 
ence 31j are shown plotted in figure 7 and compared with calculated values of 

continuum heat transfer, free -molecule 
heat transfer, and the bridged heat 
transfer value obtained from equa- 
tion ( 20 ) . The data shown are 
stagnation -point he at -transfer measure- 
ments on a spherical body at nominal 
Mach numbers of 5-7 and 8 with stagna- 
tion temperatures of 2100 ° to 2300 ° R. 
The quantities, q QC (eq. ( 15 )), q F M 
(eq. (17)), and q OQ (eq. (20)) are all 
normalized using the continuum heat- 
transfer rate with zero vorticity, 
q oc . This normalizing factor is 
vort=0 

obtained from equation (15) with the 
The Reynolds number used as abscissa. 



Figure 7.- Comparison of calculated convective 
heat transfer bridging with experiment for 
the case without blowing. 


vorticity correction dropped (Cq = 0 ) . 

Re s , is based on the enthalpy velocity (eq. (13b)) and the stagnation gas 
properties, as originally used in reference 31- 


Re s — 


%., m 

10 Ho 


(C 20 ) 


In the experiments of reference 31 * the Reynolds number was varied by mainly 
varying the free -stream density, D, and the nose radius, R. It is seen from 
the figure that the calculated bridging checks well with experiment. 

Comparisons were also made (not shown here) with measured data from 
reference 32 for subsonic heat transfer from spheres. Equation ( 20 ) checks 
these data reasonably well. 

It is concluded from the comparisons made that the bridging relation in 
equation (l8) should be a useful approximation for the general case with blow- 
ing since equation (20) is simply a specialization (with the form unchanged) 
of equation (l8), and equation (l8) has the correct asymptotes for the general 
case. Experimental verification for the blowing case would be desirable. 


Surface Shear 

For the calculation of the (x derivative of) surface shear in the 
transitional regime, we bridge between equations (35) and (36) to obtain equa- 
tion (37)' The first collision model used is the same as that used for con- 
vective heat transfer, except that we are now concerned with the transfer of 
x momentum rather than energy. 
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I 


0 



Sketch (e) 


We use the coordinate system shown 
in sketch (e), and we introduce 


(du/dx) x=o 

(du e /dx ) x=0 


u 



(C21) 


We assume that collisions occurring 
within an effective collision thick- 
ness, A t (a kind of boundary layer 
thickness), affect surface shear, while 
those occurring outside A r have a 
negligible effect. Making use of equa- 
tion (ll), we rewrite equation ( 36 ) as 


<FM = 



(C22) 


where K 4 contains the con-version of units and the accommodation coefficient, 
A cm , considered to he approximately unity. Using the same reasoning and 
analogous development as that used for convective heat transfer, we can write 


q- 1 =1 

1 W 


(sr - § 


all 5X 


rr* I — q- I 

T w t wFM 


^ 8 X(l - u*) 


all SX 

We evaluate 6 X according to equation (CIO) and obtain: 

^Aq- 


T* = T* f (1 - u*)e" Y / A S + 6 - A t/ A 

T W - t wFM / u )e A 

L u o 

We use (see sketch (e)) 


and have 


q- t — q“ ^ 

w wFM 


L w o 


u * - f w 


,&T - (B ^ f + ^ 


(C23a) 


(C23b) 


(C24) 


(C25) 


(C26) 
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As was done with the convective heat transfer bridging, we evaluate 
f(Y/A r ) in the simplest way (see sketch (e)): 



The integration of equation (C26) then yields 


»w 


=* T. 


wFM 




A 1 



(C27) 


(028) 


The ratio, A/A,., becomes very large in the free -molecule limit and very small 
in the continuum limit in a manner similar to A/A^ (as described above in the 
Convective Heat Transfer section). When A/A^. is large, r^J- in equa- 
tion ( 028 ) approaches as it should. When A/A^- is small, mast 

approach t^- c . From equation (028) we get 


A _ T wc 
a t t wFM 


(C29) 


Since the characteristic thickness, At, has not been specified, it can now be 
defined as having values that satisfy equation (C29) • We are thus forcing our 
bridging relationship to have the correct continuum asymptote. 


We can show directly that equation (C29) is a consistent relationship by 
making use of equations (5.6, 6.23, 6.25, 6.29, 6. 30, 6.32, and 6.4l) of 
reference 30. We can also show the consistency of equation (C29) by using the 
heat -transfer relations that we have. By combining equations (l6) , (17), (35), 
and (36) we obtain 


Then, using equation 


wc 


wFM 


a 3 a, 


cq 


E a 

1 + -f v 


A-cm P21 


(C19b), we have 


wc 


wFM 


Eq 

A-3&CQ ^ A + p — ^ 

l cm *J p 2i 




(C30) 


(C31) 


This gives equation (C29) when 
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A t = 


"cm 


*J P21 


A 3 A< 


•cq. 


E3 

1 + -f t 




(C32) 


When we substitute equation (C29) into (C 28 ) , we have the bridging 
relation for the x gradient of surface shear. 


T * = t * 

1 W' WC 


1 - e 


WJ 


Fm/ t -k 


WC 


( 37 ) 


If functional forms for f(Y/A T ) other than the simple one in equation (C2?) 
are chosen for the derivation given above , the bridging equations will be more 
complicated than equation (37) but similarly behaved. It is possible that 
some of the more complicated forms may furnish more accurate bridgings for 
some specific shear situations. 

Direct experimental verifications of equation ( 37 ) applied to the 
stagnation region with blowing have not been found. However, it is believed 
that this relationship is a reasonable and useful approximation for surface 
shear in the transitional regime. The expression is a monotonic function of 
t wFm/ t ^c an b i-b bas the correct free -molecule and continuum regime asymptotes. 

Although equation (37) "was derived for a stagnation region boundary 
layer, it is of possible interest to compare calculated results with subsonic 
flat plate shear measurements since these data are available. This comparison 
is shown in figure 8. The measured data are from figure 9(a) of reference 32. 

The agreement shown in figure 8 is 
thought to be surprisingly good. 

An attempt was also made to 
compare calculated shear from equa- 
tion (37) with measured shear in low 
speed Couette flow as reported in 
reference 33 (not shown) ; the agreement 
obtained was fair. An equation derived 
for Couette flow by the authors of 
reference 33 checks the data very 
closely; this illustrates that one 
should generally prefer a bridging 
equation derived for a specific situa- 
tion if this be available. 



Figure 8.- Comparison of calculated skin friction 
bridging with measured skin friction for 
subsonic flat plate (without blowing). 


Drag Bridging 

For flight calculations, we use the trajectory equations (52) to (54), 
in which the quantity, M/CpA, appears. The evaluations of Cp and M/CpA are 
given in equations (57) and (59) .> respectively, equation (57) being the 
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bridging formula for Cp. A bridging relation for the drag coefficient is 
clearly necessary for entry flights since Cp must initially have a free- 
molecule evaluation, and, for many entries, the final value of Cp will be 
the continuum value . 

To evaluate the drag bridging, we again make use of our first -collision 
model (in a treatment that amounts to a further approximation with the model). 
We again use an unspecified characteristic effective collision thickness, Ap, 
and assume that collisions outside of Ap have a negligible effect on drag. 

We assume, then, that the molecules that have a longer free path than Ap 
make their first effective collision with the body and contribute to free- 
molecule drag. The other molecules collide with each other within the depth, 
a d ; these molecules bathe the body in a continuum fluid and contribute to con- 
tinuum drag. A more rigorous development would require Ap to vary with posi- 
tion on the body, but we will take Ap to be some average value for the whole 
body. Similarly, the value of the mean free path of the molecules near the 
body will depend on position on the body, but we will use a nominal or 
averaged mean free path, A^, for the gas between the shock wave and the body. 

According to our model, we can sum the free -molecule and continuum drags 
and obtain 


C D^oo A - C DFM3oaFM A + C DC ( 3 0 oC A 


(C33a) 


J D 



+ C DC 



(C33"b) 


The dynamic pressure ratios are evaluated as density or 

SxFM _ P qqFM -Ad/A 2 
^ PqO 


ioC _ _ -A D /A 2 

Pco Poo 


mass fraction ratios: 
(C34a) 

(C34b) 


Combining equations (C 3 I 4 .) with (033^)., we have 

C D = C DC + ( C DFM " C Dc) e Al) ^ 2 


(C35) 


We can insert the quantity, E s , from equation ( 58 ) into equation (C 35 ) to 
obtain 


C P = C DC 


(l + E 9 e _A D/ A2 ^ 


(C36) 
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We observe that equations (C 35 ) and (C3 6 ) have the correct asymptotes. The 
variable, Aq/A^, is a reciprocal Knudsen number. To this point, the deriva- 
tion is similar to that found in reference 34 , but from a somewhat different 
point of view. In reference 34, Ap is specified as the shock standoff dis- 
tance, evaluated empirically in terms of p and R for a sphere; this recip- 
rocal Knudsen number is then multiplied by a constant factor to fit 
experimental data. 

We can postulate a functional form for Ap as 

= Re °°> "body shape) 

From dimensional considerations, we write 


Ap - Rfp (Re^, body shape) 


where the nose radius, R, characterizes the body size. Then we have 


A = IL r 

A 2 A 2 P 


(C37a) 


We can write 


^ = RP 21 (d 10~ 6 )ffl 
A^ ^ 2^2 


(C37^) 


This latter form is equivalent to the exponent in equation ( 7 ) of reference 34 
when P 21 ^p is taken as a constant for a spherical body shape. As shown in 
reference 34, with the constant properly adjusted, a good fit is obtained with 
the drag data for a sphere at high speeds in air and in helium. 


For practical calculations, it seems preferable to avoid evaluating A^* 
We can use the classical relation between viscosity and mean free path 
(eq. ( 119 ) of ref. 27 ) which states that viscosity is proportional to the 
product of Ap and a mean molecular speed. We can express A 2 p 2 in terms of 
an undissociated value, A^p^; this involves using fp,(pt 2 , T s ), the ratio of 
actual viscosity to the undissociated (Sutherland) viscosity. (Values of f^ 
for air are tabulated in table VI of ref. 35 .) We make use of the Sutherland 
formula. 



(C38) 


to relate to free-stream conditions, and we obtain the approximation 


^2^2 \rA 


/m 2 


co H oo \j it^q 


Tr + b(T R AbJ 

T s + b 




(C39) 
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Using equation (106) of reference 27, we can write 


. = 

AdoPcQ Q 


and put equation (C39) in the form 

KPi = 


_ /njj; 

• 2-2 Soo 7 m M 


Ta + hCTs/Tpo)' 


T s + h 


Substituting this expression into equation (C37h) we have 

( RD)lO~ 6 P ai fg 


We let 


and we have 


Aa K A m oo /“a 


T s + h(T s /Tj 


T s + h 


f 7 = 


10 s P 21 fp 


K V\x> 


Ts + ^(Ts/Too) 


moo 


T s + b 


Ad 

7\pi 


= RDf - 


(C40) 


(C4l) 


(C42) 


(C43) 


(C44) 


where 


f 7 = fy (Reoo^ Moo, gas, body shape) (C45) 

For given flight conditions and gas, f y will he a function of body shape 
(includes angle of attack). For many flight cases, body shape does not change 
much in the passage through the transitional regime; it can be expected that 
generally f y will not undergo a large variation in the transitional passage. 
We now approximate f y as a constant for the flight of a given body and use 

f 7 = 15(1 + E 14 ) (c46) 


where E X 4 can be assigned the value zero for a sphere in air, and the con- 
stant, 15 , has been selected to match experimental data as described below. 
We combine equations (C3 6 ), (C44), and (C46) to obtain our drag bridging 
equation. 


C D “ C DC 


1 -f Ec$e 


-is(Hj)(i+E 14 ) 


(57) 
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The parameter, Ei 4 , can he adjusted, in principle, to account approximately 
for the effects of all the variables in equation (C45) (during passage through 
the transitional regime). Unless data are at hand, E 14 for air would nor- 
mally be assigned the value zero, as for a sphere. Although the (presumably 
modest) dependence of E i4 on body shape may not be known, the values of Cp 
calculated from equation (57) will also depend directly on body shape through 
Cpp and Eg. The constant, 15, in equations (C46) and (57) was determined by 
matching calculated values of drag with the measured drag data for spheres in 
air of reference 3 6. The comparison is shown in figure 9* Corresponding with 

calculated values of Cp for figure 9 , 
the values of RD used in equa- 
tion (57) were converted to Re^, 
because the data of reference 36 are 
plotted with Reoo as abscissa. The 
data of reference 3 6 were obtained in 
air at a nominal settling chamber tem- 
perature of 300° K over a Mach number 
range from 3-8 to 4-3* The data indi- 
cate that any Mach number effect on 
the transition of the drag coefficient 
is probably small. The value, E 14 = 0 
was also used to obtain a good fit with 
sphere drag data in air reported in 
reference 34 (not shown) . One set of 
these data was obtained in undissociated air at a nominal settling chamber 
temperature of 25 00° K over a Mach number range from 15-96 to 20. 90. Another 
set of data was obtained in dissociated air over a range of hypersonic Mach 
numbers (11.34 to 58.7) at a nominal settling chamber temperature of 9000° K. 

A third set of sphere drag data was obtained in helium with hypersonic flow 
conditions. The helium data require E 14 ^ 2 for a good fit. 

The drag coefficient, Cp, is used only in equation (59) for the quantity 
M/CpA. Equation (59) contains an empirical evaluation of the variation of 
body mass, M, with the stagnation point surface recession, X. Any change in 
the continuum drag coefficient, Cpc, can also be accounted for in the 
empiricism of the M/CpA evaluation. 



Figure 9. - Comparison of calculated and measured 
values of the drag coefficient for spheres in 
the transitional regime. 


Example of Free -Molecule -Continuum Bridging 

The bridging relations developed in this appendix are essentially 
approximations to be used for the conditions for which they were derived. 

For other conditions, more rigorously developed equations may be available, 
as, for example, for cases of shear in Couette flow such as reported in refer- 
ence 33* Under conditions for which no well developed bridging relations 
exist, the formulas given in this appendix are recommended as engineering 
approximations. The structure of these equations insures that the equations 
have the correct asymptotes, and this tends to limit inaccuracies. 
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Representative curves obtained 
from the bridging equations are pre- 
sented in figure 10. The ordinate 
quantities, normalized to their free- 
molecule values, are plotted against a 
reciprocal Knudsen number which was 
varied by varying the nose radius, R. 
Although the bridging equations are 
general, the numerical inputs to the 
equations, and therefore the curves 
obtained, will depend on the individual 
case calculated. The calculated curves 
in the figure are for a tektite glass. 

A curve for the quantity q^ w (eq. (l8)) 
is not shown; this quantity is the 
product of the quantities, \| r and q QQ 
(see eq. (22b)), for which the curves are shown. The dotted lines on the 
right side of the figure are the continuum regime asymptotes for the various 
quantities, while on the left side the free -molecule asymptote for all quanti- 
ties is unity. It is noted that the various quantities approach their asymp- 
totes at different rates. (Under some conditions, one can consider that all of 
the quantities plotted are not even in the same regime.) The figure illus- 
trates probably the most important feature of the bridging relations used. 

They automatically place control of the various quantities in the appropriate 
controlling regime, the free -molecule, transitional, or continuum. 

Curves corresponding to those of figure 10 can be plotted for other 
materials, and most of the curves will be similar, but somewhat displaced. 

The heat transfer curve (c^ ) is approximately universal; the drag curve (Cp) 
is influenced considerably by body shape; the other curves (t^., v w , \|r) will 
vary somewhat with relative rates of vaporization (or reaction) of the 
material being considered. 



.01 .1 1.0 10 100 1000 
BODY RADIUS — , R/ 

MEAN FREE PATH 

Figure 10.- Representative calculation of 
bridging between free -molecule and 
continuum flow. 
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APPENDIX D 


USE OF COMPUTING PROGRAM 


The computing program can be used to solve a variety of problems 
involving surface-type ablation, fata are arbitrarily read in subject to the 
limitations given below in the listing of Input Data. As explained in the 
section METHOD OF THE NUMERICAL PROGRAM, the stability parameter, Z, should be 
< l/2 for all grid points and at all times. The quantity, Z, is printed out 
by the program so that its values may be observed. 


Computing Program Options 

The numerical computing program has six major groups of options as listed 
below. 

1. Running conditions 

(a) Normal wind tunnel, KF = 1. 

(b) Rarefied wind tunnel, KF = 2. 

(includes all wind tunnel cases, but computing time is longer than 
with option (a)) 

(c) Flight, KF = 3* 

Both of the wind-tunnel options, (a) and (b), allow wind-tunnel condi- 
tions to be changed once, if desired, and also for the wind tunnel to be shut 
off. At time, t c ^ a (GAMAI), the free-stream density, enthalpy velocity, and 
free-stream velocity are changed, respectively, to De^a (C 7 ), V c ] ia (C 8 ), and 
Voocha (OKBAR). At time, t Q ff (OMCO), the wind tunnel is shut off, and the 
calculations are continued while the model cools. If these changes are not 
desired, the values of t c ;k a and t 0 ^ can be set larger than the time corre- 
sponding to the final time line number, NE (see Time Sketch below). 

2. Internal radiation 

(a) Transparent, KG = 1. 

(b) Opaque, KG = 2. 

(c) Semitransparent, KG = 3* 

3 . Surface conditions 

(a) Evaporation or sublimation, KCH = 1. 

(b) Surface chemical reaction, KCH = 2. 


62 



4. Initial conditions for flight case (KF = 3) > (see appendix B) . 

(a) Normal entry into an atmosphere using a computed exponential tempera- 

ture profile in the body. With an assumed T wi (> T 0 ), the initial 
values, Dj_ and Aq are computed by the program, KDEL = 1. 

(b) Arbitrary initial values of atmospheric density, Dq, and thermal 

thickness of exponential temperature profile in the body, Aq, 

KDEL = 2. 

5. Back face boundary conditions 

(a) Back face aerodynamically exposed, KBAK = 1. 

(b) Backing material forming a heat sink, KBAK = 2. (For a heat sink of 

zero heat capacity, or an adiabatic back boundary, KBAK = 1 should 
be used. ) 

6. Planet and atmosphere for flight case (KF = 3) 

(a) Earth entry with the ARDC atmosphere (approximated exponentially with 

3 programmed values of scale height) . Earth radius programmed at 
6440 km (R p in eq. (53)), KC5 = 1. 

(b) Arbitrary planet with exponential atmosphere having arbitrary scale 

height (initial, two intermediate, and final values), KC5 = 2. 


Nomenclature of Computing Program 

The nomenclature used in the computing program is in symbolic FORTRAN 
language. Separate listings of input and output data are shown below. 


Input Data 

Input data are listed below in their order of card punching. Actual 
card formats are shown in the Input Card Format Sketch. All input data are 
printed out by the program in an initial readout (see Sample Case, below). 

A quantity listed as an option is defined in the section. Computing Program 
Options, above. 

Following the definition of a quantity, the value of an option selection 
may be shown. The particular quantity is not needed (and not used) for other 
values of the option selection. Unused quantities are normally assigned the 
value zero; in any case, the input card formats must be maintained. 

The maximum number of grid points to be used in the finite difference 
spacing is 98 ; this is the maximum value for the quantity, MF (see Spacing 
Sketch) . 
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W 


CARO 


4 

6 

12 

16 

20 24 

26 

32 

56 

.0 

44 

40 

52 

56 

60 

_KF 

KG 

Nl 

NZ 

NF KZ 

K3 

K4 

LN 

M2 

M3 

MF 

ji 

J2 

KC5 

M3 

KMI 

KCM3 

It* 3 

(CKf| (III 

M2 

(CHI 

(M2 

ttlf 

(CH 

(DEL 

(6(( 




ALL NUMBERS IN 14 FORMAT (MUST BE RIGHT JUSTIFIED) 



Input Card Format Sketch 


Ml, M2, M3 ARE BREAK POINTS FOR CHANGING A 77 


FINE SPACED 
e.g.,LN=4 


M.HIII 




(A^ ) 


(A??,) 




GRID POINT NO.=M I 


e.g.,K2=4 (CALCULATED) 

e.g.,Mi = i7 


MF 

(MAX = 98) 


Nl, N2 ARE BREAK POINTS FOR CHANGING At 


| At , = DDT I 

TIME LINE NQ = N 1 
t = tj 


- A t g = DDT 2 - 


H 


At jl At 2 , At 3 , ARE INCREMENTS USED IN FINITE DIFFERENCE EQUATION 

Time Sketch 


ONE TIME LINE 

KMI, KM2 , KM3 ARE BREAK POINTS FOR PRINTING INCREMENTS 


(KCM2) 


(KCM3) 


GRID POINT NO. I KMI KM2 KM3 MF 

KCMI, KCM2, KCM3, KCMF ARE GRID-POINT PRINTING INCREMENTS 


INITIALLY AY, = Ar) i 

Atj 2 , Ai 7 3 , At ; 4 ARE INCREMENTS USED IN FINITE DIFFERENCE EQUATION 
A 17 , INCREMENT USED FOR INTERPOLATIONS 


BETWEEN TIME LINES 

KNI, KN2 ARE BREAK POINTS FOR PRINTING INCREMENTS 


(KCN I) 


(KCN2) 


(KCNF) 


INITIAL VALUE OF LENGTH/A 7?, < 1665 REQUIRED FOR 
TRANSPARENT CASE (KG= I) 


TIME LINE NO. I KNI KN2 NF 

KCNI, KCN 2, KCNF ARE TIME LINE PRINTING INCREMENTS 


Spacing Sketch 


Printing Sketch 


Card A 

(All numbers are integers in 1 4 FORMAT) 

KF Running condition option 

KG Internal radiation option 

Nl Time line number at which the finite time increment At (DTT) changes 
from Ati (DDTl) to At 2 (DDT2) (see Time Sketch). 

N2 Time line number at which the finite time increment At (DTT) changes 
from At 2 (DDT2) to At 3 (DDT3) (see Time Sketch). 

NF Final time line number (see Time Sketch). 

K2 Defined by Arj 2 ~ (K2) At^ (or Ay 2 = (K2) Ay 1 ) . Can be 1 for opaque 

and semitransparent cases (KG =2, 3), but must be at least 2 and even 
for transparent case (KG = l) (see Spacing Sketch). 
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K3 

K4 

LN 

M2 

M3 

MF 1 

J1 

J2 

KC5 

KCT 

KM1 

KM2 

KM3 

KCMl 

KCM3 

KN1 

KN2 


Defined by Ar| 3 = (K3) (or Ay 3 = (K3) Ay 2 ) . Can be 1 (see 

Spacing Sketch) . 

Defined by At ] 4 = (K4) Arj 3 (or Ay 4 = (K4) Ay 3 ) . Can be 1 (see 
Spacing Sketch) . 

Increments of At ] 2 spacing over which Ar^ spacing exists. 

Must be > 4 (see Spacing Sketch). 

Grid point at which space increment changes from Arj 2 to At ] 3 
(see Spacing Sketch). 

Grid point at which space increment changes from Ar| 3 to Ar| 4 
(see Spacing Sketch). 

Grid point at back face. Maximum value = 98 (see Spacing Sketch). 

Order of interpolation for T (TE) in the fine spaced (At] 1 ) 
region. 

Order of interpolation for A^ (YDEL) . 

Planet and atmosphere option for flight case (KF = 3) . 

Maximum number of iterations to determine front face temperature, 
T v , within allowable error selected (ALLOW) . If exceeded, 
calculation will stop. 

Grid point at which printing interval on one time line changes 
from KCMl to KCM2 (see Printing Sketch). 

Grid point at which printing interval on one time line changes 
from KCM2 to KCM3 (see Printing Sketch). 


Card B 

(All numbers are integers in 1 4 FORMAT) 

Grid point at which printing interval on one time line changes 
from KCM3 to KCMF (see Printing Sketch). 

j- Printing intervals of grid points (see Printing Sketch). 

Time line number at which time line printing interval changes 
from KCN1 to KCN2 (see Printing Sketch). 

Time line number at which time line printing interval changes 
from KCN2 to KCNF (see Printing Sketch). 
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Time line printing intervals (see Printing Sketch). 


KCN 1 ,KCN 2 , 

KCNF 

KCH Surface condition option. 


KDEL 

Initial condition option for flight case (KF = 3) . 


KBAK 

Back face boundary condition option. 


For cards 1 - 11 , there are normally 8 numbers per card, each number 
E9.3 FORMAT. 



Card 1 



A 1 

Ai, equation (27) 

B 1 

P 


A 2 

A 2 , equations (lO),(ll) 

B 2 

B 2 , equation ( 84 ) 


A 3 

A3 , equations (32), (3 5) 

B 3 

B3, equation ( 84 ) 


A 4 

A4 , equation (15) 

B 4 

B 4 , equation (85) 




Card 2 



B 5 

B 5 , equation (85) 

B 9 

B 9 , equation (87) 


b6 

B e , equation (86) 

BIO 

M, equation ( 29 ) 


B 7 

B T , equation (86) 

Bll 

Bn, equations (28), 

(30) 

B8 

T 

- L o 

B 12 





Card 3 



B 13 

epp for opaque case (KG = 

2 ) > € max 

for semitransparent case 

(KG 


Bl 4 Bi 4 , equation (85) 

B 15 e-gj, for opaque and semitransparent cases (KG = 2, 3) 

Bl 6 B 1S , equations ( 46 ), ( 47 ), ( 51 ); B i6 = 1 for transparent case (KG = l); 
B 1S = 0 for opaque and semitransparent cases (KG = 2 , 3) 

Cl Ci, equation (60) 

C 2 C 2 , equation (60) 

C 3 C 3 , equations (5 6), (59) 

C 4 C 4 , equations (56), ( 59 ) 
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I 



C5 

C6 


Card 4 


Sh x ^ for flight case (KF = 3), and (KC5 = 2). 

Cq, equation (15) 

C7 L/D r , equation (53) for flight case (KF = 3) '> for wind tunnel 

cases (KF = 1, 2). 

C8 gp/10 5 , equation ( 53"b ) for flight case (KF = 3); V cila for wind 

tunnel cases (KF =1., 2) . 

El Ei, equation (86) 

E2 E 2 , equation (87) 

E3 E 3 , equation (88) 

E4 E 4 , equation (44-) 


Card 5 


E5 

E 5 , equation (44) 

E9 

Eg , equation (58) 

E 6 

E 6 , equation (44) 

E10 

E10, equation (88) 

E7 

Ey, equation (26) 

Ell 

Eh, equation (6l) 

E8 

E 8 , equations (34), (35) 

E12 

E12, equation (6l) 



Card 6 


E13 

E 13 , equation (50) 



El4 

E14, equations (57), (59) 



E15 

Eis, equation (lOa) 



E16 

Eie, equation (12) 




SIGMA or = 1.369X10" 12 

RO R^, equation (60) 

OMCO (M/GqqA)^, equation (59) for flight case (KF = 3); ^off f° r ' w ’ in( ^ 
tunnel cases (KF = 1, 2) . 

OKBAR K, equation (12) for flight case (KF = 3)) V , for wind tunnel 
cases (KF = 1, 2) . 
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Card 7 


CHI1 y, x , equation (12) 

VCINFI V ooi for flight case (KF = 3); for wind tunnel cases (KF = 1, 2), 

V OT = VCIHFI until t > t cha or t > t off . 

GAMAI 7j_ for flight case (KF = 3); t cha for wind tunnel cases (KF = 1, 2). 

TWI T w: jj > T 0 for entry flight case (KF - 3> KDEL = l) . 

DELTW Initial guess for incrementing T wi in iterating for Tw in the 
second time line (can he zero) . 

ALL0W e aZ 

HS h s for wind tunnel cases (KF = 1, 2) ; constant at read -in value 

until t > t c k a or t > t Q ff . 

DC I Dj_ for flight case (KF =3); read in only when KDEL = 2; with 

KDEL = 1, Dj_ is calculated by the program. For wind tunnel cases 
(KF = ly 2), read in; D = DC I until t > t c ^ a or t > t 0 ff. 

Card 8 

DDT1 Ati (see Time Sketch). 

DDT2 At 2 (see Time Sketch). 

LDT3 At 3 ( see Time Sketch) . 

DY1 Atj 1 = initial value of A (see Spacing Sketch); initial value of 

length/AT ) 1 < 1665 is required for the transparent case (KG = l) . 

RH021 P21; read in for wind tunnel cases (KF =1, 2) ; for flight case 

(KF = 3), P21 is calculated continuously (eq. (6 l)). 

0M0 %/F' e l ua ii° n (48); used only when back face is aerodynamically 

exposed (KBAK = l) . 

EI7 E 17 , equation (60) 

El8 E 1s , equations (56)>(59)> used for flight case (KF = 3). 

Card 9 

DELTJ Initial guess for incrementing T w in iteration for time lines 3, 

4, 5; should not be zero. 

E19 equation (49) ] used when hack face is aerodynamically exposed 

(KBAK = 1) . 
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RB 

R-b, equations (48), ( 50 ); used when hack face is aero dynamic ally exposed 
(KBAK = 1) . 

E20 

E 2Q , equation ( 87 ) 


E21 

A C q, equation (l7)> (KF = 2, 3 ); if unknown use A C q = 1. 


E22 

A crr i ; equation ( 36 ), (KF = 2, 3 ); if unknown use A cm = 1. 


E23 

Ac V , equations (39c, 39d) , (KCH = l) ] if unknown, evaluate using 
A cv = 1 in equation (39c) . 

E24 

B, equation (39e), (KCH = 2). 



Card 10 


E25 

Efji, equation (39c), (KCH =■ 2). 


E 26 

c, equation (51), (KBAK = 2); (if c = 0, use KBAK = 1 and 0M0 = 0) . 

E27 

Rp, equation (53), flight case (KF =3), and (KC5 = 2 ). 


E28 

j? , flight case (KF =3), and (KC5 = 2 ). 


E29 

pT^ , flight case (KF = 3), and (KC5 = 2 ). 


E30 

Sb 2 , flight case (KF = 3), and (KC5 = 2 ). 


E 31 

S^g, flight case (KF = 3 ), and (KC5 = 2 ). 


E32 

£>h±> in equation (BIO), flight case (KF =3); used only when KDEL = 1 
and KC 5 = 2 ; has arbitrary value, but may equal Sj ll (C5) • 


Card 11 


E33 

E 33 , equation (45) , semitransparent case (KG = 3)* 


E 34 

Aj_, in equation (Bl) ; read in only for flight case (KF = 3)j with 
arbitrary initial conditions (KDEL = 2) . Otherwise the program 
computes Aq from equation (B2b) . 

E 35 

E 35 , equation ( 28 ) 


E 36 

E 36 , transparent case (KG =1); equation (97); E 3e = 1.0 has 
good energy balances. 

given 

E37 

E 37 , transparent case (KG = l) ; equation ( 98 ); E 37 =• 1.0 has 
good energy balances. 

given 

E 38 

E39,E40 

E 38 , flight case (KF = 3); equation ( 55 ). 
(Open) 
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Card 12 


This card is only required for the transparent case (KG = l) . The format 
for this card is shown in the Input Card Format Sketch. The FORTRAN quanti- 
ties ALPHA, ALPHA2, EMAX, RN are in FORMAT E9-3; IT is in FORMAT 14. 

ALPHA <x, equation (4l) 

ALPHA2 a 2 , equation (4l) 

EMAX e maX j equation (42) 

RN n, equations (41, 42) 

IT Spacing of calculation and printing of F and g = bF/bj in the 

fine spaced region (see Spacing Sketch). IT must he an integer 
factor of K2 hut must he less than K2; for example, if K2 = 6, 

IT can he 1, 2, or 3* IT = 3 means that F and g are printed 
out for every third grid point in the fine spaced region; in the 
remaining regions they are printed for every grid point. The 
value of IT affects somewhat the accuracy of the finite differ- 
ence and energy balance calculations. IT = 1 gives the greatest 
accuracy. A larger value of IT decreases accuracy while reducing 
computing time. 


Output Data 

All input data are printed out in an initial output format (see Sample 

Case helow) . In addition, the following quantities (some calculated) are 

printed in an initial output (see Sample Case helow). 

Q01I q vi 

Q02I q^j_ 

QOI ^wRi 

TTI t^; equation (b4) for wind tunnel cases (KF = 1, 2) ; equals zero 

for flight case (KF = 3) • 

DELI A^; equation (B2h) for wind tunnel cases (KF = 1, 2) and flight 

case (KF = 3) with KDEL = 1; for flight case with KDEL = 2, Aj[ 
equals E34. 

DCI Dj_; read in quantity for wind tunnel cases (KF = 1, 2) and flight 

case (KF =3) with KDEL = 2; from equation (BIO) for flight case 
with KDEL = 1. 

RHO = D/1226 (for this printing D = Dj_) . 

Ml Ml = 1 + (K2)(LN); see Spacing Sketch. 
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j! , 

( The spacing of printing intervals is arbitrarily selected (see subsection 
Input Data and the Printing Sketch) . For each time line printed out, the 
following quantities appear. (See section. Sample Case, below for format.) 

N N, time line number 

I 

| TT t 

I 

I R E 

? SNGMA sin 7 


RHO 

Pco 

PDP 

P" 

VC 

V 

P 

P 

VCINF 

V 

00 

PHI 

Cp 

PSI 

V 

A2D 

AaD 

PSIBAR 


FBAR 

F 

QOC 

^OC 

FINT 

F t 

QOFM 

q FM 

ADGABS 

a Dg 

QO 

% 

AGABS 

a g 

Q00 

^oo 

EPSRER 

€ bf 

PSQOO 

q tw 

EPSEN 

€ ff 

CK3U 


FRS 

f rs 

PSQOF 

q v 

PT2 

Pts 

QR 


COIN 

^con 

TDEL 

A 

RIN 

Sr ad 

YDEL 


VAP 

^vap 

TMDL 

). 

STOR 

^stor 

RH021 

P 21 

VC ON 

Q 

ycon 

TAUCP 

T* 

WC 

UCON 

^ucon 

TAUFMP 

T wFM 

RESID 

Sres 

TAUWP 

T f 

1 V 

ERR 

e 

rr 


OMC M/CpA. 

XIP v sr = dX/dt 

CHI X 

DTT At 
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DTDYW 

(Ml 

VARG 

^vcon 


YTNT 

^-ucon 

PSQO 

q 

DRES 

^res 

con 



FRW 

^rad 

KC 

KC, number 


required 

DVAP 

Q 



^-vap 

TW 

Tv 

DSTOR 

‘Is tor' exce P t at t = H> 

WC 

I*; I 


when DSTOR = (stored 

1 v vc 1 


energy at tjJ/At; stored 

VWFM 

Iv 1 


energy at t± = absolute 



value of second term in 

VW 

, 


equation ( 82 d) . 

v w 


iterations 
obtain T w . 


Also for each time line, the following arrayed quantities print out in 
columns for the selected print spacing of grid points (see Printing Sketch 
and Sample Case below for format). 

M M (grid point number) 

Y y 

TE T 

UDB u 

VBAR v 

XRAT X r 

0 g 

FR F 

Z Z, must be < 1/2 for stability. 

Sample Case 

To illustrate the use of the computing program, a sample case has been 
selected which describes the entry of a transparent tektite into the Earth's 
atmosphere. Figure 11 shows the input data for the sample case; figure 12 
shows the printing out of the input data and the initial calculated values; 
figure 13 shows the printing out for a typical time line. 
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Figure 11.- Input data 
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Figure I 3 .- Print-out of a typical time line for sample case 
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